diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
index a403098..a4818f0 100644
--- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
+++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
@@ -1,373 +1,378 @@
// Copyright 2016-2021 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.SetXTitle("cos#theta_{#mu}-p_{#mu} (GeV)");
+ fSettings.SetYTitle("d^{2}#sigma/dp_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)");
fSettings.SetAllowedTypes("FULL,DIAG/FREE,SHAPE,FIX/SYSTCOV/STATCOV",
"FIX/FULL");
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();
+ fSaveFine = false;
};
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;
fMode = event->Mode;
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() {
// Do standard conversion.
Measurement1D::ConvertEventRates();
// 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 (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++;
}
}
if (fMCHist_Modes) {
fMCHist_Modes->Reset();
for (std::map >::iterator it =
fMCModeHists_Slices.begin();
it != fMCModeHists_Slices.end(); ++it) {
bincount = 0;
for (size_t i = 0; i < nangbins; i++) {
it->second[i]->Scale(fScaleFactor / (angular_binning_costheta[i + 1] -
angular_binning_costheta[i]),
"width");
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
fMCHist_Modes->SetBinContent(it->first, bincount + 1, 0, 0,
it->second[i]->GetBinContent(j + 1));
fMCHist_Modes->SetBinError(it->first, bincount + 1, 0, 0,
it->second[i]->GetBinError(j + 1));
bincount++;
}
}
}
}
return;
}
void T2K_CC0pi_XSec_2DPcos_nu_I::FillMCSlice(double x, double y, double w) {
if (fMCHist_Modes && !fMCModeHists_Slices.count(fMode)) {
for (size_t i = 0; i < nangbins; ++i) {
std::stringstream ss("");
ss << "T2K_CC0pi_XSec_2DPcos_nu_I_MODE_ " << fMode << "_slice" << i
<< std::endl;
fMCModeHists_Slices[fMode].push_back(
static_cast(fMCHist_Slices[i]->Clone(ss.str().c_str())));
fMCModeHists_Slices[fMode].back()->Reset();
}
}
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);
if (fMCHist_Modes) {
fMCModeHists_Slices[fMode][i]->Fill(x, w);
}
}
}
}
void T2K_CC0pi_XSec_2DPcos_nu_I::SetHistograms() {
// Read in 1D Data Histograms
TFile input(
(FitPar::GetDataBase() + "/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root")
- .c_str(),
- "READ");
+ .c_str(), "READ");
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);
+ "T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D",
+ 400, 0.0, 30.0, 100, -1.0, 1.0);
+ fMCHist_Fine2D->GetXaxis()->SetTitle("p_{#mu} (GeV)");
+ fMCHist_Fine2D->GetYaxis()->SetTitle("cos#theta_{#mu}");
+ fMCHist_Fine2D->GetZaxis()->SetTitle(fSettings.GetYTitle().c_str());
+
fMCHist_Fine2D->SetDirectory(NULL);
SetAutoProcessTH1(fMCHist_Fine2D);
TH2D *tempcov = (TH2D *)input.Get("analysis1_totcov");
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);
}
}
// Read in 2D Data Slices and Make MC Slices
int bincount = 0;
for (size_t i = 0; i < nangbins; i++) {
fDataHist_Slices.push_back(
(TH1D *)input.Get(Form("dataslice_%lu", i))->Clone());
fDataHist_Slices[i]->SetNameTitle(
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);
+ fDataHist_Slices.back()->GetXaxis()->SetTitle("p_{#mu} (GeV)");
+ fDataHist_Slices.back()->GetYaxis()->SetTitle(fSettings.GetYTitle().c_str());
// 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++;
}
}
assert(bincount == tempcov->GetNbinsX());
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)));
}
}
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%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);
- fDataHist =
- new TH1D("T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D",
- "T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D", bincount, 0, bincount);
+ fDataHist = new TH1D("T2K_CC0pi_XSec_2DPcos_nu_I_data_1D",
+ ("T2K_CC0pi_XSec_2DPcos_nu_I_data_1D"+fSettings.PlotTitles()).c_str(),
+ 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);
}
SetShapeCovar();
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%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%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();
}
}
}
diff --git a/src/T2K/T2K_CC0pi_XSec_H2O_2DPcos_anu.cxx b/src/T2K/T2K_CC0pi_XSec_H2O_2DPcos_anu.cxx
index bb7dfd5..5a3ba37 100755
--- a/src/T2K/T2K_CC0pi_XSec_H2O_2DPcos_anu.cxx
+++ b/src/T2K/T2K_CC0pi_XSec_H2O_2DPcos_anu.cxx
@@ -1,190 +1,199 @@
// Copyright 2016-2021 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_H2O_2DPcos_anu.h"
static size_t nmombins = 7;
static double mom_binning[] = { 400., 530., 670., 800., 1000., 1380., 2010., 3410. };
static int ncosbins[] = { 2, 3, 3, 3, 3, 3, 2 };
static double costheta_binning[][7] = { { 0.84, 0.94, 1. },
{ 0.85, 0.92, 0.96, 1. },
{ 0.88, 0.93, 0.97, 1. },
{ 0.90, 0.94, 0.97, 1. },
{ 0.91, 0.95, 0.97, 1. },
{ 0.92, 0.96, 0.98, 1. },
{ 0.95, 0.98, 1. } };
//********************************************************************
T2K_CC0pi_XSec_H2O_2DPcos_anu::T2K_CC0pi_XSec_H2O_2DPcos_anu(nuiskey samplekey) {
//********************************************************************
+ // Don't use the FINE histogram
+ fSaveFine = false;
+
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CC0pi_XSec_H2O_2DPcos_anu. \n"
"Target: H2O \n"
"Flux: T2K 2.5 degree off-axis (ND280-P0D) \n"
"Signal: CC0pi\n"
"arXiv:1908.10249";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
- fSettings.SetXTitle("cos#theta_{#mu}-p_{#mu}");
+ fSettings.SetXTitle("cos#theta_{#mu}-p_{#mu} (GeV)");
fSettings.SetYTitle("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("H,O");
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pi_XSec_H2O_2DPcos_anu");
fSettings.DefineAllowedSpecies("numub");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/water molecule
fScaleFactor = ((GetEventHistogram()->Integral("weight")/(fNEvents+0.)) * 18E-38 / (TotalIntegratedFlux("")));
// Setup Histograms
SetHistograms();
// Final setup
FinaliseMeasurement();
};
bool T2K_CC0pi_XSec_H2O_2DPcos_anu::isSignal(FitEvent *event){
return SignalDef::isT2K_CC0piAnuP0D(event, EnuMin, EnuMax);
};
void T2K_CC0pi_XSec_H2O_2DPcos_anu::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();
double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect()));
fXVar = pmu;
fYVar = CosThetaMu;
return;
};
void T2K_CC0pi_XSec_H2O_2DPcos_anu::FillHistograms(){
Measurement1D::FillHistograms();
if (Signal){
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_H2O_2DPcos_anu::ConvertEventRates(){
for (size_t i = 0; i < nmombins; i++){
fMCHist_Slices[i]->GetSumw2();
}
// Do standard conversion.
Measurement1D::ConvertEventRates();
for (size_t i = 0; i < nmombins; ++i) {
fMCHist_Slices[i]->Scale(1. / (mom_binning[i + 1] - mom_binning[i]));
}
// Now Convert into 1D list
fMCHist->Reset();
int bincount = 0;
for (size_t i = 0; i < nmombins; i++){
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++){
fMCHist->SetBinContent(bincount+1, fMCHist_Slices[i]->GetBinContent(j+1));
bincount++;
}
}
return;
}
void T2K_CC0pi_XSec_H2O_2DPcos_anu::FillMCSlice(double x, double y, double w){
for (size_t i = 0; i < nmombins; ++i) {
if ((x > mom_binning[i]) && (x <= mom_binning[i + 1])) {
fMCHist_Slices[i]->Fill(y, w);
}
}
}
void T2K_CC0pi_XSec_H2O_2DPcos_anu::SetHistograms(){
// Read in 1D Data Histograms
fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/AntiNuMuH2O/AntiNuMuOnH2O_unreg.root").c_str(),"READ");
+ fPlotTitles = fSettings.GetFullTitles();
// Read in 1D Data
fDataHist = (TH1D*) fInputFile->Get("xsecDataRelease");
fDataHist->SetNameTitle((fName + "_data").c_str(),
(fName + "_data" + fPlotTitles).c_str());
int Nbins = fDataHist->GetNbinsX();
// Read relative covariance matrix
TH2D* tempcov = (TH2D*) fInputFile->Get("covDataRelease");
// Make absolute covariance matrix
fFullCovar = new TMatrixDSym(Nbins);
for (int i = 0; i < Nbins; i++){
for (int j = 0; j < Nbins; j++){
(*fFullCovar)(i,j) = tempcov->GetBinContent(i+1, j+1)*fDataHist->GetBinContent(i+1)*fDataHist->GetBinContent(j+1)*1E76;
}
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
// Read in 2D Data Slices and Make MC Slices
int bincount = 0;
for (uint i = 0; i < nmombins; ++i) {
- fDataHist_Slices.push_back(new TH1D(Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_data_Slice%u",i),Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_data_Slice%u",i),ncosbins[i],costheta_binning[i]));
+ fDataHist_Slices.push_back(new TH1D(Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_data_Slice%u",i),
+ Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_data_Slice%u",i),
+ ncosbins[i],costheta_binning[i]));
+ fDataHist_Slices[i]->GetXaxis()->SetTitle("cos#theta_{#mu}");
+ fDataHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetYTitle().c_str());
+
for (int j = 0; j < ncosbins[i]; ++j) {
fDataHist->SetBinError(bincount+1,sqrt((*fFullCovar)(bincount,bincount))*1E-38);
fDataHist_Slices[i]->SetBinContent(j+1, fDataHist->GetBinContent(bincount+1));
fDataHist_Slices[i]->SetBinError(j+1, sqrt((*fFullCovar)(bincount,bincount))*1E-38);
bincount++;
}
fMCHist_Slices.push_back((TH1D*) fDataHist_Slices[i]->Clone());
fMCHist_Slices[i]->SetNameTitle(Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_MC_Slice%u",i),
- (Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_MC_Slice%u",i)));
+ (Form("T2K_CC0pi_XSec_H2O_2DPcos_anu_MC_Slice%u",i)));
SetAutoProcessTH1(fDataHist_Slices[i],kCMD_Write);
SetAutoProcessTH1(fMCHist_Slices[i]);
}
return;
};
diff --git a/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx b/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx
index 06fbefb..96f4f52 100644
--- a/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx
+++ b/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx
@@ -1,257 +1,254 @@
#include "T2K_CCinc_XSec_2DPcos_nu_nonuniform.h"
#include "T2K_SignalDef.h"
// ***********************************
// Implemented by Alfonso Garcia, Barcelona (now NIKHEF)
// Clarence Wret, Rochester
// (Alfonso was the T2K analyser)
// ***********************************
//********************************************************************
T2K_CCinc_XSec_2DPcos_nu_nonuniform::T2K_CCinc_XSec_2DPcos_nu_nonuniform(
nuiskey samplekey) {
//********************************************************************
+ fSaveFine = false;
fAllowedTypes += "/GENIE/NEUT";
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CCinc_XSec_2DPcos_nu_nonuniform sample. \n"
"Target: CH \n"
"Flux: T2K FHC numu \n"
"Signal: CC-inclusive \n"
"https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.012004";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
- fSettings.SetXTitle("");
- // fSettings.SetXTitle("p_{#mu} (GeV)");
- // fSettings.SetYTitle("cos#theta_{#mu}");
+ fSettings.SetXTitle("p_{#mu}-cos#theta_{#mu}");
fSettings.SetYTitle("#frac{d^{2}#sigma}{dp_{#mu}dcos#theta_{#mu}} "
- "[#frac{cm^{2}}{nucleon/GeV/c}]");
+ "(cm^{2}/nucleon/GeV)");
fSettings.SetEnuRange(0.0, 30.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K CC-inclusive p_{#mu} cos#theta_{#mu}");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * 1E-38 / fNEvents /
TotalIntegratedFlux();
// Default to using the NEUT unfolded data
UnfoldWithGENIE = false;
// Check option
if (fSettings.Found("type", "GENIE"))
UnfoldWithGENIE = true;
// Tell user what's happening
if (UnfoldWithGENIE) {
NUIS_LOG(SAM, fName << " is using GENIE unfolded data. Want NEUT? Specify "
"type=\"NEUT\" in your config file");
} else {
NUIS_LOG(SAM, fName << " is using NEUT unfolded data. Want GENIE? Specify "
"type=\"GENIE\" in your config file");
}
// Setup Histograms
SetHistograms();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
// Signal is simply a CC inclusive without any angular/momentum cuts
bool T2K_CCinc_XSec_2DPcos_nu_nonuniform::isSignal(FitEvent *event) {
return SignalDef::isCCINC(event, 14, EnuMin, EnuMax);
};
void T2K_CCinc_XSec_2DPcos_nu_nonuniform::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;
};
// Fill up the MCSlice
void T2K_CCinc_XSec_2DPcos_nu_nonuniform::FillHistograms() {
if (Signal)
FillMCSlice(fXVar, fYVar, Weight);
}
// Modification is needed after the full reconfigure to move bins around
void T2K_CCinc_XSec_2DPcos_nu_nonuniform::ConvertEventRates() {
// Do standard conversion.
Measurement1D::ConvertEventRates();
// First scale MC slices also by their width in Y
fMCHist_Slices[0]->Scale(1.0 / 0.25);
fMCHist_Slices[1]->Scale(1.0 / 0.50);
fMCHist_Slices[2]->Scale(1.0 / 0.20);
fMCHist_Slices[3]->Scale(1.0 / 0.15);
fMCHist_Slices[4]->Scale(1.0 / 0.11);
fMCHist_Slices[5]->Scale(1.0 / 0.09);
fMCHist_Slices[6]->Scale(1.0 / 0.07);
fMCHist_Slices[7]->Scale(1.0 / 0.05);
fMCHist_Slices[8]->Scale(1.0 / 0.04);
fMCHist_Slices[9]->Scale(1.0 / 0.025);
fMCHist_Slices[10]->Scale(1.0 / 0.015);
// Now Convert into 1D list
fMCHist->Reset();
int bincount = 0;
for (int i = 0; i < nSlices; i++) {
for (int j = 0; j < fMCHist_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++;
}
}
};
void T2K_CCinc_XSec_2DPcos_nu_nonuniform::FillMCSlice(double x, double y,
double w) {
if (y >= -1.0 && y < -0.25)
fMCHist_Slices[0]->Fill(x, w);
else if (y >= -0.25 && y < 0.25)
fMCHist_Slices[1]->Fill(x, w);
else if (y >= 0.25 && y < 0.45)
fMCHist_Slices[2]->Fill(x, w);
else if (y >= 0.45 && y < 0.6)
fMCHist_Slices[3]->Fill(x, w);
else if (y >= 0.6 && y < 0.71)
fMCHist_Slices[4]->Fill(x, w);
else if (y >= 0.71 && y < 0.80)
fMCHist_Slices[5]->Fill(x, w);
else if (y >= 0.80 && y < 0.87)
fMCHist_Slices[6]->Fill(x, w);
else if (y >= 0.87 && y < 0.92)
fMCHist_Slices[7]->Fill(x, w);
else if (y >= 0.92 && y < 0.96)
fMCHist_Slices[8]->Fill(x, w);
else if (y >= 0.96 && y <= 0.985)
fMCHist_Slices[9]->Fill(x, w);
else if (y >= 0.985 && y <= 1.0)
fMCHist_Slices[10]->Fill(x, w);
};
void T2K_CCinc_XSec_2DPcos_nu_nonuniform::SetHistograms() {
// Read in 1D Data Histograms
- TFile *fInputFile =
- new TFile((FitPar::GetDataBase() +
+ TFile *fInputFile = new TFile((FitPar::GetDataBase() +
"T2K/CCinc/nd280data-numu-cc-inc-xs-on-c-2018/histograms.root")
- .c_str(),
- "OPEN");
+ .c_str(), "OPEN");
// Number of theta slices in the release
nSlices = 11;
// Data release includes unfolding with NEUT or GENIE as prior
// Choose whichever the user specifies
std::string basename;
if (UnfoldWithGENIE)
basename = "hist_xsec_data_prior_neut_cthbin";
else
basename = "hist_xsec_data_prior_neut_cthbin";
// Read in 2D Data Slices and Make MC Slices
// Count the number of bins we have in total so we can match covariance matrix
int bincount = 0;
for (int i = 0; i < nSlices; i++) {
// Get Data Histogram
- // fDataHist_Slices.push_back((TH1D*)fInputFile->Get(Form("dataslice_%i",i))->Clone());
fDataHist_Slices.push_back(
(TH1D *)fInputFile->Get(Form("%s%i", basename.c_str(), i))
->Clone(
- Form("T2K_CCinc_XSec_2DPcos_nu_nonuniform_slice%i_data", i)));
+ Form("T2K_CCinc_XSec_2DPcos_nu_nonuniform_Slice%i_data", i)));
fDataHist_Slices[i]->SetDirectory(0);
fDataHist_Slices[i]->Scale(1E-39);
- fDataHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str());
+ fDataHist_Slices[i]->GetXaxis()->SetTitle("p_{#mu} (GeV)");
+ fDataHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetYTitle().c_str());
// Count up the bins
bincount += fDataHist_Slices.back()->GetNbinsX();
// Make MC Clones
fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone(
Form("T2K_CCinc_XSec_2DPcos_nu_nonuniform_Slice%i_MC", i)));
fMCHist_Slices[i]->Reset();
fMCHist_Slices[i]->SetDirectory(0);
fMCHist_Slices[i]->SetLineColor(kRed);
- fMCHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str());
+ fMCHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetYTitle().c_str());
SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write);
SetAutoProcessTH1(fMCHist_Slices[i]);
}
fDataHist =
new TH1D((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str(), bincount, 0, bincount);
fDataHist->SetDirectory(0);
int counter = 0;
for (int i = 0; i < nSlices; ++i) {
// Set a nice title
std::string costitle = fDataHist_Slices[i]->GetTitle();
costitle = costitle.substr(costitle.find("-> ") + 3, costitle.size());
std::string found = costitle.substr(0, costitle.find(" < "));
std::string comp = costitle.substr(costitle.find(found) + found.size() + 3,
costitle.size());
comp = comp.substr(comp.find(" < ") + 3, comp.size());
costitle = "cos#theta_{#mu}=" + found + "-" + comp;
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
fDataHist->SetBinContent(counter + 1,
fDataHist_Slices[i]->GetBinContent(j + 1));
fDataHist->SetBinError(counter + 1,
fDataHist_Slices[i]->GetBinError(j + 1));
// Set a nice axis
if (j == 0)
fDataHist->GetXaxis()->SetBinLabel(
counter + 1, Form("%s, p_{#mu}=%.1f-%.1f", costitle.c_str(),
fDataHist_Slices[i]->GetBinLowEdge(j + 1),
fDataHist_Slices[i]->GetBinLowEdge(j + 2)));
else
fDataHist->GetXaxis()->SetBinLabel(
counter + 1,
Form("p_{#mu}=%.1f-%.1f", fDataHist_Slices[i]->GetBinLowEdge(j + 1),
fDataHist_Slices[i]->GetBinLowEdge(j + 2)));
counter++;
}
}
// The correlation matrix
// Major in angular bins, minor in momentum bins: runs theta1, pmu1, pmu2,
// theta2, pmu1, pmu2, theta3, pmu1, pmu2 etc The correlation matrix
TH2D *tempcov = NULL;
if (UnfoldWithGENIE)
tempcov = (TH2D *)fInputFile->Get("covariance_matrix_genie");
else
tempcov = (TH2D *)fInputFile->Get("covariance_matrix_neut");
fFullCovar = new TMatrixDSym(bincount);
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
(*fFullCovar) *= 1E38 * 1E38;
(*covar) *= 1E-38 * 1E-38;
fInputFile->Close();
};