diff --git a/data/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root b/data/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root
index d61ec36..d1d4028 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 e80c5ba..d43d0c6 100644
--- a/data/T2K/CC0pi/makedatafile_t2kcc0pi.py
+++ b/data/T2K/CC0pi/makedatafile_t2kcc0pi.py
@@ -1,296 +1,297 @@
from ROOT import *
from array import *
+from math 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")
datapoly = TH2Poly("datapoly","datapoly", 0.0,30.0, -1.0, 1.0)
hist = None
binedges = []
histedgeslist = []
xsecvals = []
histxseclist = []
binlimits = [3,8,15,22,30,39,47,58,67]
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] ) + 1
-
+
xval = round(float(GetLowEdge( data[2] )),4)
yval = round(float(GetLowEdge( data[1] )),4)
xhig = round(float(GetHighEdge( data[2] )),4)
yhig = round(float(GetHighEdge( data[1] )),4)
-
+
xsec = float( data[3] ) * 1E-38
datapoly.AddBin( xval, yval, xhig, yhig )
datapoly.SetBinContent( datapoly.GetNumberOfBins(), xsec)
binedges.append( xval )
xsecvals.append( xsec )
- if ibin in binlimits:
+ if ibin in binlimits:
binedges.append( xhig )
histedgeslist.append(binedges)
histxseclist.append(xsecvals)
binedges = []
xsecvals = []
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)
# Get Covariances (keep in 1E-38 cm^2) \
nbins = 67
statcov = TH2D("analysis1_statcov","analysis1_statcov", nbins, 0.0, float(nbins), nbins, 0.0, float(nbins));
systcov = TH2D("analysis1_systcov","analysis1_systcov", nbins, 0.0, float(nbins), nbins, 0.0, float(nbins));
normcov = TH2D("analysis1_normcov","analysis1_normcov", nbins, 0.0, float(nbins), nbins, 0.0, float(nbins));
totcov = TH2D("analysis1_totcov","analysis1_totcov", nbins, 0.0, float(nbins), nbins, 0.0, float(nbins));
with open("covariance_statisticUncertainty_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
xi, yi = GetIndex(data[0])
cov = float(data[1])
statcov.SetBinContent(xi + 1, yi + 1, cov)
with open("covariance_shapeSystematics_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
xi, yi = GetIndex(data[0])
cov = float(data[1])
systcov.SetBinContent(xi + 1, yi + 1, cov)
with open("covariance_fluxNormalizationSystematics_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
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)
data1D = TH1D("datahist","datahist", datapoly.GetNumberOfBins(), 0.0, float(datapoly.GetNumberOfBins()));
for i in range(datapoly.GetNumberOfBins()):
data1D.SetBinContent(i+1, datapoly.GetBinContent(i+1));
data1D.SetBinError(i+1, sqrt(totcov.GetBinContent(i+1,i+1))*1E-38)
outfile.cd()
for i, obj in enumerate(histedgeslist):
print obj
hist = TH1D("dataslice_" + str(i), "dataslice_" + str(i), len(obj)-1, array('f',obj))
for j in range(hist.GetNbinsX()):
hist.SetBinContent(j+1, histxseclist[i][j])
- hist.GetXaxis().SetRangeUser(obj[0], obj[len(obj)-2])
+ hist.GetXaxis().SetRangeUser(obj[0], obj[len(obj)-1])
hist.Draw("HIST")
gPad.Update()
hist.SetNameTitle("dataslice_" + str(i),"dataslice_" + str(i))
hist.Write()
outfile.cd()
datahist.Write()
counthist.Write()
maphist.Write()
datapoly.Write()
data1D.Write()
statcov.Write()
systcov.Write()
totcov.Write()
normcov.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] ) + 1
xval = GetMiddle( data[2] )
yval = GetMiddle( data[1] )
xsec = float( data[3] ) * 1E-38
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 + 1, yi + 1, 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()
+totcov.Write()
+normcov.Write()
# LocalWords: xval
diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
index ea0706d..882cc20 100644
--- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
+++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
@@ -1,270 +1,328 @@
// 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 "T2K_SignalDef.h"
#include "T2K_CC0pi_XSec_2DPcos_nu_I.h"
+static size_t nangbins = 9;
+static double angular_binning_costheta[] = {-1, 0, 0.6, 0.7, 0.8,
+ 0.85, 0.9, 0.94, 0.98, 1};
+
//********************************************************************
T2K_CC0pi_XSec_2DPcos_nu_I::T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CC0pi_XSec_2DPcos_nu_I sample. \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0pi\n"
"https://journals.aps.org/prd/abstract/10.1103/"
"PhysRevD.93.112012 Analysis I";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("P_{#mu} (GeV)");
fSettings.SetYTitle("cos#theta_{#mu}");
fSettings.SetZTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)");
fSettings.SetAllowedTypes("FULL,DIAG/FREE,SHAPE,FIX/SYSTCOV/STATCOV",
"FIX/FULL");
- fSettings.SetEnuRange(0.0, 10.0);
fSettings.DefineAllowedTargets("C,H");
+ fSettings.SetEnuRangeFromFlux(fFluxHist);
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pi_XSec_2DPcos_nu_I");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
+ fMaskMomOverflow = false;
+ if (samplekey.Has("mask_mom_overflow")) {
+ fMaskMomOverflow = samplekey.GetB("mask_mom_overflow");
+ }
+
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) *
1E-38 / (TotalIntegratedFlux()));
+ assert(std::isnormal(fScaleFactor));
+
// Setup Histograms
SetHistograms();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
bool T2K_CC0pi_XSec_2DPcos_nu_I::isSignal(FitEvent *event) {
return SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, SignalDef::kAnalysis_I);
};
void T2K_CC0pi_XSec_2DPcos_nu_I::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double pmu = Pmu.Vect().Mag() / 1000.;
double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect()));
fXVar = pmu;
fYVar = CosThetaMu;
return;
};
void T2K_CC0pi_XSec_2DPcos_nu_I::FillHistograms() {
Measurement1D::FillHistograms();
if (Signal) {
fMCHist_Fine2D->Fill(fXVar, fYVar, Weight);
FillMCSlice(fXVar, fYVar, Weight);
}
}
// 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_I::ConvertEventRates() {
-
- for (int i = 0; i < 9; i++) {
- fMCHist_Slices[i]->GetSumw2();
- }
-
// Do standard conversion.
Measurement1D::ConvertEventRates();
- // First scale MC slices also by their width in Y
- fMCHist_Slices[0]->Scale(1.0 / 1.00);
- fMCHist_Slices[1]->Scale(1.0 / 0.60);
- fMCHist_Slices[2]->Scale(1.0 / 0.10);
- fMCHist_Slices[3]->Scale(1.0 / 0.10);
- fMCHist_Slices[4]->Scale(1.0 / 0.05);
- fMCHist_Slices[5]->Scale(1.0 / 0.05);
- fMCHist_Slices[6]->Scale(1.0 / 0.04);
- fMCHist_Slices[7]->Scale(1.0 / 0.04);
- fMCHist_Slices[8]->Scale(1.0 / 0.02);
+ // Scale MC slices by their bin area
+ for (size_t i = 0; i < nangbins; ++i) {
+ fMCHist_Slices[i]->Scale(fScaleFactor / (angular_binning_costheta[i + 1] -
+ angular_binning_costheta[i]),
+ "width");
+ }
// Now Convert into 1D list
fMCHist->Reset();
int bincount = 0;
- for (int i = 0; i < 9; i++) {
+ for (size_t i = 0; i < nangbins; i++) {
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
fMCHist->SetBinContent(bincount + 1,
fMCHist_Slices[i]->GetBinContent(j + 1));
fMCHist->SetBinError(bincount + 1, fMCHist_Slices[i]->GetBinError(j + 1));
bincount++;
}
}
return;
}
void T2K_CC0pi_XSec_2DPcos_nu_I::FillMCSlice(double x, double y, double w) {
- if (y >= -1.0 and y < 0.0)
- fMCHist_Slices[0]->Fill(x, w);
- else if (y >= 0.0 and y < 0.6)
- fMCHist_Slices[1]->Fill(x, w);
- else if (y >= 0.6 and y < 0.7)
- fMCHist_Slices[2]->Fill(x, w);
- else if (y >= 0.7 and y < 0.8)
- fMCHist_Slices[3]->Fill(x, w);
- else if (y >= 0.8 and y < 0.85)
- fMCHist_Slices[4]->Fill(x, w);
- else if (y >= 0.85 and y < 0.90)
- fMCHist_Slices[5]->Fill(x, w);
- else if (y >= 0.90 and y < 0.94)
- fMCHist_Slices[6]->Fill(x, w);
- else if (y >= 0.94 and y < 0.98)
- fMCHist_Slices[7]->Fill(x, w);
- else if (y >= 0.98 and y <= 1.00)
- fMCHist_Slices[8]->Fill(x, w);
+ for (size_t i = 0; i < nangbins; ++i) {
+ if ((y >= angular_binning_costheta[i]) &&
+ (y < angular_binning_costheta[i + 1])) {
+ fMCHist_Slices[i]->Fill(x, w);
+ }
+ }
}
void T2K_CC0pi_XSec_2DPcos_nu_I::SetHistograms() {
// Read in 1D Data Histograms
- fInputFile = new TFile(
+ TFile input(
(FitPar::GetDataBase() + "/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root")
.c_str(),
"READ");
- // fInputFile->ls();
-
- // Read in 1D Data
- fDataHist = (TH1D *)fInputFile->Get("datahist");
-
fMCHist_Fine2D = new TH2D("T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D",
"T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D", 400, 0.0, 30.0,
100, -1.0, 1.0);
+ fMCHist_Fine2D->SetDirectory(NULL);
SetAutoProcessTH1(fMCHist_Fine2D);
- TH2D *tempcov = (TH2D *)fInputFile->Get("analysis1_totcov");
+ TH2D *tempcov = (TH2D *)input.Get("analysis1_totcov");
- fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX());
- for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
- for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
+ fFullCovar = new TMatrixDSym(tempcov->GetNbinsX());
+ for (int i = 0; i < tempcov->GetNbinsX(); i++) {
+ for (int j = 0; j < tempcov->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
- covar = StatUtils::GetInvert(fFullCovar);
- fDecomp = StatUtils::GetDecomp(fFullCovar);
-
- // Read in 2D Data
- fDataPoly = (TH2Poly *)fInputFile->Get("datapoly");
- fDataPoly->SetNameTitle("T2K_CC0pi_XSec_2DPcos_nu_I_datapoly",
- "T2K_CC0pi_XSec_2DPcos_nu_I_datapoly");
- SetAutoProcessTH1(fDataPoly, kCMD_Write);
- fDataHist->Reset();
// Read in 2D Data Slices and Make MC Slices
int bincount = 0;
- for (int i = 0; i < 9; i++) {
+ for (size_t i = 0; i < nangbins; i++) {
- // Get Data Histogram
- // fInputFile->ls();
fDataHist_Slices.push_back(
- (TH1D *)fInputFile->Get(Form("dataslice_%i", i))->Clone());
+ (TH1D *)input.Get(Form("dataslice_%lu", i))->Clone());
+
fDataHist_Slices[i]->SetNameTitle(
- Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%i", i),
- (Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%i", i)));
+ Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%lu", i),
+ (Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%lu, cos(#theta) [%f,%f] ",
+ i, angular_binning_costheta[i],
+ angular_binning_costheta[i + 1])));
+ fDataHist_Slices.back()->SetDirectory(NULL);
// Loop over nbins and set errors from covar
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
fDataHist_Slices[i]->SetBinError(
j + 1, sqrt((*fFullCovar)(bincount, bincount)) * 1E-38);
+ bincount++;
+ }
+ }
- // std::cout << "Setting data hist " <<
- // fDataHist_Slices[i]->GetBinContent(j+1) << " " <<
- // fDataHist_Slices[i]->GetBinError(j+1) << std::endl;
- fDataHist->SetBinContent(bincount + 1,
- fDataHist_Slices[i]->GetBinContent(j + 1));
- fDataHist->SetBinError(bincount + 1,
- fDataHist_Slices[i]->GetBinError(j + 1));
+ assert(bincount == tempcov->GetNbinsX());
- bincount++;
+ if (fMaskMomOverflow) {
+ MaskMomOverflow();
+ bincount = fFullCovar->GetNcols();
+ }
+
+ std::vector> data_slice_bcbes;
+ for (size_t i = 0; i < nangbins; i++) {
+ for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
+ data_slice_bcbes.push_back(
+ std::make_pair(fDataHist_Slices[i]->GetBinContent(j + 1),
+ fDataHist_Slices[i]->GetBinError(j + 1)));
}
+ }
- // Make MC Clones
+ for (size_t i = 0; i < nangbins; i++) {
fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone());
+ fMCHist_Slices.back()->SetDirectory(NULL);
+ fMCHist_Slices.back()->Reset();
+ fMCHist_Slices.back()->SetLineColor(kRed);
fMCHist_Slices[i]->SetNameTitle(
- Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%i", i),
- (Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%i", i)));
+ Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%lu", i),
+ (Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%lu, cos(#theta) [%f,%f] ", i,
+ angular_binning_costheta[i], angular_binning_costheta[i + 1])));
+ }
+
+ covar = StatUtils::GetInvert(fFullCovar);
+ fDecomp = StatUtils::GetDecomp(fFullCovar);
- SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write);
- SetAutoProcessTH1(fMCHist_Slices[i]);
- // fMCHist_Slices[i]->Reset();
+ fDataHist =
+ new TH1D("T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D",
+ "T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D", bincount, 0, bincount);
+ fDataHist->SetDirectory(NULL);
+ for (size_t i = 0; i < data_slice_bcbes.size(); ++i) {
+ fDataHist->SetBinContent(i + 1, data_slice_bcbes[i].first);
+ fDataHist->SetBinError(i + 1, data_slice_bcbes[i].second);
}
+ fMCHist = (TH1D *)fDataHist->Clone("T2K_CC0pi_XSec_2DPcos_nu_I_MC_1D");
+ fMCHist->SetDirectory(NULL);
return;
-};
+}
+
+void T2K_CC0pi_XSec_2DPcos_nu_I::MaskMomOverflow() {
+
+ std::vector bins_to_cut;
+ size_t nallbins = 0;
+ for (size_t i = 0; i < nangbins; i++) {
+
+ std::vector slice_bin_edges;
+ slice_bin_edges.push_back(
+ fDataHist_Slices[i]->GetXaxis()->GetBinLowEdge(1));
+ for (int j = 0; j < (fDataHist_Slices[i]->GetNbinsX() - 1); j++) {
+ slice_bin_edges.push_back(
+ fDataHist_Slices[i]->GetXaxis()->GetBinUpEdge(j + 1));
+ nallbins++;
+ }
+
+ bins_to_cut.push_back(nallbins++);
+ TH1D *tmp = new TH1D(fDataHist_Slices[i]->GetName(),
+ fDataHist_Slices[i]->GetTitle(),
+ slice_bin_edges.size() - 1, slice_bin_edges.data());
+ tmp->SetDirectory(NULL);
+ for (int j = 0; j < (fDataHist_Slices[i]->GetNbinsX() - 1); j++) {
+ tmp->SetBinContent(j + 1, fDataHist_Slices[i]->GetBinContent(j + 1));
+ tmp->SetBinError(j + 1, fDataHist_Slices[i]->GetBinError(j + 1));
+ }
+
+ delete fDataHist_Slices[i];
+ fDataHist_Slices[i] = tmp;
+ }
+
+ TMatrixDSym *tmpcovar = new TMatrixDSym(nallbins - bins_to_cut.size());
+ int icut = 0;
+ for (int ifull = 0; ifull < fFullCovar->GetNcols(); ifull++) {
+ if (std::find(bins_to_cut.begin(), bins_to_cut.end(), ifull) !=
+ bins_to_cut.end()) {
+ continue;
+ }
+ int jcut = 0;
+ for (int jfull = 0; jfull < fFullCovar->GetNcols(); jfull++) {
+ if (std::find(bins_to_cut.begin(), bins_to_cut.end(), jfull) !=
+ bins_to_cut.end()) {
+ continue;
+ }
+ (*tmpcovar)[icut][jcut] = (*fFullCovar)[ifull][jfull];
+ jcut++;
+ }
+ icut++;
+ }
+ delete fFullCovar;
+ fFullCovar = tmpcovar;
+}
void T2K_CC0pi_XSec_2DPcos_nu_I::Write(std::string drawopt) {
this->Measurement1D::Write(drawopt);
+ for (size_t i = 0; i < nangbins; i++) {
+ fMCHist_Slices[i]->Write();
+ fDataHist_Slices[i]->Write();
+ }
+
if (fResidualHist) {
std::vector MCResidual_Slices;
size_t tb_it = 0;
for (size_t i = 0; i < fMCHist_Slices.size(); ++i) {
- std::string name = Form("T2K_CC0pi_XSec_2DPcos_nu_I_RESIDUAL_Slice%i", i);
+ std::string name =
+ Form("T2K_CC0pi_XSec_2DPcos_nu_I_RESIDUAL_Slice%lu", i);
MCResidual_Slices.push_back(
static_cast(fMCHist_Slices[i]->Clone(name.c_str())));
MCResidual_Slices.back()->Reset();
for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) {
double bc = fResidualHist->GetBinContent(tb_it + 1);
MCResidual_Slices.back()->SetBinContent(j + 1, bc);
tb_it++;
}
MCResidual_Slices.back()->Write();
}
}
if (fChi2LessBinHist) {
std::vector MCChi2LessBin_Slices;
size_t tb_it = 0;
for (size_t i = 0; i < fMCHist_Slices.size(); ++i) {
std::string name =
- Form("T2K_CC0pi_XSec_2DPcos_nu_I_Chi2NMinusOne_Slice%i", i);
+ Form("T2K_CC0pi_XSec_2DPcos_nu_I_Chi2NMinusOne_Slice%lu", i);
MCChi2LessBin_Slices.push_back(
static_cast(fMCHist_Slices[i]->Clone(name.c_str())));
MCChi2LessBin_Slices.back()->Reset();
for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) {
double bc = fChi2LessBinHist->GetBinContent(tb_it + 1);
MCChi2LessBin_Slices.back()->SetBinContent(j + 1, bc);
tb_it++;
}
MCChi2LessBin_Slices.back()->Write();
}
}
-}
\ No newline at end of file
+}
diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h
index a1f63b6..2570f1e 100644
--- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h
+++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h
@@ -1,74 +1,72 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef T2K_CC0PI_2DPCOS_NU_NONUNIFORM_H_SEEN
#define T2K_CC0PI_2DPCOS_NU_NONUNIFORM_H_SEEN
#include "Measurement1D.h"
#include "TH2Poly.h"
#include "MeasurementVariableBox2D.h"
class T2K_CC0pi_XSec_2DPcos_nu_I : public Measurement1D {
public:
T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey);
/// Virtual Destructor
~T2K_CC0pi_XSec_2DPcos_nu_I() {};
/// Numu CC0PI Signal Definition
///
- /// /item
+ /// /item
bool isSignal(FitEvent *nvect);
/// Read histograms in a special way because format is different.
/// Read from FitPar::GetDataBase()+"/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root"
void SetHistograms();
/// Bin Tmu CosThetaMu
void FillEventVariables(FitEvent* customEvent);
// Fill Histograms
void FillHistograms();
/// Have to do a weird event scaling for analysis 1
void ConvertEventRates();
void Write(std::string drawopt);
/// \brief Create Q2 Box to save correction info
inline MeasurementVariableBox* CreateBox(){ return new MeasurementVariableBox2D(); };
private:
bool fIsSystCov, fIsStatCov, fIsNormCov;
+ bool fMaskMomOverflow;
- TH2Poly* fDataPoly;
- TH2Poly* fMCPoly;
-
- TFile* fInputFile;
TH2D* fMCHist_Fine2D;
std::vector fMCHist_Slices;
std::vector fDataHist_Slices;
void FillMCSlice(double x, double y, double w);
-
+ void MaskMomOverflow();
+
};
-
+
#endif