Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F10881359
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
32 KB
Subscribers
None
View Options
diff --git a/src/T2K/T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos.cxx b/src/T2K/T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos.cxx
index 09baa19..c5b030c 100644
--- a/src/T2K/T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos.cxx
+++ b/src/T2K/T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos.cxx
@@ -1,250 +1,249 @@
// 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_SignalDef.h"
#include "T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos.h"
static size_t nangbins = 9;
static double angular_binning_costheta[] = {-1, 0.2, 0.6, 0.7, 0.8,
0.85, 0.9, 0.94, 0.98, 1 };
//********************************************************************
T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos::T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos(nuiskey samplekey) {
//********************************************************************
// Samples overview ---------------------------------------------------
fSettings = LoadSampleSettings(samplekey);
std::string name = fSettings.GetS("name");
std::string descrip = "";
// This has to deal with NuMu FHC, and AntiNuMu RHC
if (!name.compare("T2K_NuMu_CC0pi_CH_XSec_2DPcos")){
descrip = name +". \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0pi\n"
- "arXiv:2002.09323";
+ "DOI:10.1103/PhysRevD.101.112001";
fSettings.SetTitle("T2K_NuMu_CC0pi_CH_XSec_2DPcos");
fSettings.DefineAllowedSpecies("numu");
NuPDG = 14;
LepPDG = 13;
}
else if (!name.compare("T2K_AntiNuMu_CC0pi_CH_XSec_2DPcos")){
descrip = name +". \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0pi\n"
- "arXiv:2002.09323";
-
+ "DOI:10.1103/PhysRevD.101.112001";
fSettings.SetTitle("T2K_AntiNuMu_CC0pi_CH_XSec_2DPcos");
fSettings.DefineAllowedSpecies("numub");
NuPDG = -14;
LepPDG = -13;
}
// Setup common settings
fSettings.SetDescription(descrip);
fSettings.SetXTitle("p_{#mu}-cos#theta_{#mu}");
fSettings.SetYTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)");
fSettings.SetAllowedTypes("DIAG,FULL/FREE,SHAPE,FIX/SYSTCOV/STATCOV","FIX");
fSettings.SetEnuRangeFromFlux(fFluxHist);
fSettings.DefineAllowedTargets("C,H");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * 1E-38 / (TotalIntegratedFlux()));
// Setup Histograms
SetHistograms();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
bool T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos::isSignal(FitEvent *event){
return SignalDef::isCC0pi(event, NuPDG, EnuMin, EnuMax);
};
void T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos::FillEventVariables(FitEvent* event){
if (event->NumFSParticle(LepPDG) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(LepPDG)->fP;
double pmu = Pmu.Vect().Mag()/1000.;
double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect()));
fXVar = pmu;
fYVar = CosThetaMu;
return;
};
void T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos::FillHistograms(){
Measurement1D::FillHistograms();
if (Signal){
FillMCSlice( fXVar, fYVar, Weight );
}
}
void T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos::ConvertEventRates(){
for (int i = 0; i < nangbins; i++){
if(NuPDG==14) fMCHistNuMu_Slices[i]->GetSumw2();
else if(NuPDG==-14) fMCHistAntiNuMu_Slices[i]->GetSumw2();
}
// Do standard conversion.
Measurement1D::ConvertEventRates();
// Scale MC slices by their bin area
for (size_t i = 0; i < nangbins; ++i) {
if(NuPDG==14) fMCHistNuMu_Slices[i]->Scale(1. / (angular_binning_costheta[i + 1] - angular_binning_costheta[i]));
else if(NuPDG==-14) fMCHistAntiNuMu_Slices[i]->Scale(1. / (angular_binning_costheta[i + 1] - angular_binning_costheta[i]));
}
// Now Convert into 1D lists
fMCHist->Reset();
int bincount = 0;
for (int i = 0; i < nangbins; i++){
if(NuPDG==14){
for (int j = 0; j < fDataHistNuMu_Slices[i]->GetNbinsX(); j++){
fMCHist->SetBinContent(bincount+1, fMCHistNuMu_Slices[i]->GetBinContent(j+1));
bincount++;
}
}
else if(NuPDG==-14){
for (int j = 0; j < fMCHistAntiNuMu_Slices[i]->GetNbinsX(); j++){
fMCHist->SetBinContent(bincount+1, fMCHistAntiNuMu_Slices[i]->GetBinContent(j+1));
bincount++;
}
}
}
return;
}
void T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos::FillMCSlice(double x, double y, double w){
for (size_t i = 0; i < nangbins; ++i) {
if ((y > angular_binning_costheta[i]) && (y <= angular_binning_costheta[i + 1])) {
if(NuPDG==14) fMCHistNuMu_Slices[i]->Fill(x, w);
else if(NuPDG==-14) fMCHistAntiNuMu_Slices[i]->Fill(x, w);
}
}
}
void T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos::SetHistograms(){
// Read in 1D Data Histograms
fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/JointNuMu-AntiNuMu/JointNuMuAntiNuMuCC0piXsecDataRelease.root").c_str(),"READ");
TH1D* hLinearResult;
if(NuPDG==14) hLinearResult = (TH1D*) fInputFile->Get("hNuMuCC0piXsecLinearResult");
else if(NuPDG==-14) hLinearResult = (TH1D*) fInputFile->Get("hAntiNuMuCC0piXsecLinearResult");
int Nbins = hLinearResult->GetNbinsX();
std::string histoLinearNuType;
if(NuPDG==14) histoLinearNuType = "NuMuCC0pi";
else if(NuPDG==-14) histoLinearNuType = "AntiNuMuCC0pi";
// Now Convert into 1D list
fDataHist = new TH1D(("LinarResult" + histoLinearNuType).c_str(),("LinarResult" + histoLinearNuType).c_str(),Nbins,0,Nbins);
for (int bin = 0; bin < Nbins; bin++){
fDataHist->SetBinContent(bin+1, hLinearResult->GetBinContent(bin+1));
}
// Make covariance matrix
fFullCovar = new TMatrixDSym(Nbins);
TMatrixDSym* tmpcovstat = (TMatrixDSym*) fInputFile->Get("JointNuMuAntiNuMuCC0piXsecCovMatrixStat");
TMatrixDSym* tmpcovsyst = (TMatrixDSym*) fInputFile->Get("JointNuMuAntiNuMuCC0piXsecCovMatrixSyst");
for(int ibin=0; ibin<Nbins; ibin++){
for(int jbin=0; jbin<Nbins; jbin++){
if(NuPDG==14) (*fFullCovar)(ibin,jbin) = ((*tmpcovstat)(ibin,jbin) + (*tmpcovsyst)(ibin,jbin))*1E76;
else if(NuPDG==-14) (*fFullCovar)(ibin,jbin) = ((*tmpcovstat)(ibin+Nbins,jbin+Nbins) + (*tmpcovsyst)(ibin+Nbins,jbin+Nbins))*1E76;
}
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
fDataHist->Reset();
int bincount = 0;
for (int i = 0; i < nangbins; i++){
if(NuPDG==14){
// Make slices for data
fDataHistNuMu_Slices.push_back((TH1D*)fInputFile->Get(Form("hXsecNuMuCC0piDataSlice_%i",i))->Clone());
fDataHistNuMu_Slices[i]->SetNameTitle(Form("T2K_NuMu_CC0pi_2DPcos_data_Slice%i",i),
(Form("T2K_NuMu_CC0pi_2DPcos_data_Slice%i",i)));
// Loop over nbins and set errors from covar
for (int j = 0; j < fDataHistNuMu_Slices[i]->GetNbinsX(); j++){
fDataHistNuMu_Slices[i]->SetBinError(j+1, sqrt((*fFullCovar)(bincount,bincount))*1E-38);
fDataHist->SetBinContent(bincount+1, fDataHistNuMu_Slices[i]->GetBinContent(j+1));
fDataHist->SetBinError(bincount+1, fDataHistNuMu_Slices[i]->GetBinError(j+1));
bincount++;
}
// Save MC slices
fMCHistNuMu_Slices.push_back((TH1D*) fDataHistNuMu_Slices[i]->Clone());
fMCHistNuMu_Slices[i]->SetNameTitle(Form("T2K_NuMu_CC0pi_2DPcos_MC_Slice%i",i), (Form("T2K_NuMu_CC0pi_2DPcos_MC_Slice%i",i)));
SetAutoProcessTH1(fDataHistNuMu_Slices[i],kCMD_Write);
SetAutoProcessTH1(fMCHistNuMu_Slices[i]);
}
else if(NuPDG==-14) {
// Make slices for data
fDataHistAntiNuMu_Slices.push_back((TH1D*)fInputFile->Get(Form("hXsecAntiNuMuCC0piDataSlice_%i",i))->Clone());
fDataHistAntiNuMu_Slices[i]->SetNameTitle(Form("T2K_AntiNuMu_CC0pi_2DPcos_data_Slice%i",i),
(Form("T2K_AntiNuMu_CC0pi_2DPcos_data_Slice%i",i)));
//Loop over nbins and set errors from covar
for (int j = 0; j < fDataHistAntiNuMu_Slices[i]->GetNbinsX(); j++){
fDataHistAntiNuMu_Slices[i]->SetBinError(j+1, sqrt((*fFullCovar)(bincount,bincount))*1E-38);
fDataHist->SetBinContent(bincount+1, fDataHistAntiNuMu_Slices[i]->GetBinContent(j+1));
fDataHist->SetBinError(bincount+1, fDataHistAntiNuMu_Slices[i]->GetBinError(j+1));
bincount++;
}
// Save MC slices
fMCHistAntiNuMu_Slices.push_back((TH1D*) fDataHistAntiNuMu_Slices[i]->Clone());
fMCHistAntiNuMu_Slices[i]->SetNameTitle(Form("T2K_AntiNuMu_CC0pi_2DPcos_MC_Slice%i",i), (Form("T2K_AntiNuMu_CC0pi_2DPcos_MC_Slice%i",i)));
SetAutoProcessTH1(fDataHistAntiNuMu_Slices[i],kCMD_Write);
SetAutoProcessTH1(fMCHistAntiNuMu_Slices[i]);
}
}
return;
};
diff --git a/src/T2K/T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint.cxx b/src/T2K/T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint.cxx
index e6b2f71..f01e6a5 100644
--- a/src/T2K/T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint.cxx
+++ b/src/T2K/T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint.cxx
@@ -1,155 +1,155 @@
#include "T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint.h"
//********************************************************************
T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint::T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint(nuiskey samplekey){
//********************************************************************
fSettings = LoadSampleSettings(samplekey);
std::string descrip = "T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint. \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0pi\n"
- "arXiv:2002.09323";
+ "DOI:10.1103/PhysRevD.101.112001";
fSettings.SetTitle("T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint");
fSettings.DefineAllowedSpecies("numu, numub");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("p_{#mu}-cos#theta_{#mu}");
fSettings.SetYTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)");
fSettings.SetAllowedTypes("DIAG,FULL/FREE,SHAPE,FIX/SYSTCOV/STATCOV","FIX");
fSettings.SetEnuRange(0.0, 30.0);
fSettings.DefineAllowedTargets("C,H");
FinaliseSampleSettings();
if (fSubInFiles.size() != 2) {
NUIS_ABORT("T2K NuMu-AntiNuMu joint requires input files in format: NuMu and AntiNuMu");
}
std::string inFileNuMu = fSubInFiles.at(0);
std::string inFileAntiNuMu = fSubInFiles.at(1);
// Create some config keys
nuiskey NuMuKey = Config::CreateKey("sample");
NuMuKey.SetS("input", inFileNuMu);
NuMuKey.SetS("type", fSettings.GetS("type"));
NuMuKey.SetS("name", "T2K_NuMu_CC0pi_CH_XSec_2DPcos");
NuMuCC0pi = new T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos(NuMuKey);
nuiskey AntiNuMuKey = Config::CreateKey("sample");
AntiNuMuKey.SetS("input", inFileAntiNuMu);
AntiNuMuKey.SetS("type", fSettings.GetS("type"));
AntiNuMuKey.SetS("name", "T2K_AntiNuMu_CC0pi_CH_XSec_2DPcos");
AntiNuMuCC0pi = new T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos(AntiNuMuKey);
// Sort out the data hist
this->CombineDataHists();
// Set the covariance
SetCovariance();
// Add to chain for processing
fSubChain.clear();
fSubChain.push_back(NuMuCC0pi);
fSubChain.push_back(AntiNuMuCC0pi);
// This saves information from the sub-measurements
fSaveSubMeas = true;
FinaliseMeasurement();
};
//********************************************************************
void T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint::SetCovariance(){
//********************************************************************
fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/JointNuMu-AntiNuMu/JointNuMuAntiNuMuCC0piXsecDataRelease.root").c_str(),"READ");
TMatrixDSym* tmpcovstat = (TMatrixDSym*) fInputFile->Get("JointNuMuAntiNuMuCC0piXsecCovMatrixStat");
TMatrixDSym* tmpcovsyst = (TMatrixDSym*) fInputFile->Get("JointNuMuAntiNuMuCC0piXsecCovMatrixSyst");
fFullCovar = new TMatrixDSym(*tmpcovstat);
(*fFullCovar)+=(*tmpcovsyst);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
ScaleCovar(1E76);
return;
}
//********************************************************************
void T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint::CombineDataHists(){
//********************************************************************
TH1D *hNuMuData = (TH1D*)NuMuCC0pi->GetDataHistogram();
TH1D *hAntiNuMuData = (TH1D*)AntiNuMuCC0pi->GetDataHistogram();
int nbins = hNuMuData->GetNbinsX() + hAntiNuMuData->GetNbinsX();
fDataHist = new TH1D((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str(), nbins, 0, nbins);
fDataHist->SetDirectory(0);
int count = 0;
for (int x=0; x<hNuMuData->GetNbinsX(); ++x){
fDataHist->SetBinContent(count+1, hNuMuData->GetBinContent(x+1));
fDataHist->SetBinError(count+1, hNuMuData->GetBinError(x+1));
fDataHist->GetXaxis()->SetBinLabel(count+1, Form("NuMu CC0pi %.1f-%.1f", hNuMuData->GetXaxis()->GetBinLowEdge(x+1), hNuMuData->GetXaxis()->GetBinUpEdge(x+1)));
count++;
}
for (int x=0; x<hAntiNuMuData->GetNbinsX(); ++x){
fDataHist->SetBinContent(count+1, hAntiNuMuData->GetBinContent(x+1));
fDataHist->SetBinError(count+1, hAntiNuMuData->GetBinError(x+1));
fDataHist->GetXaxis()->SetBinLabel(count+1, Form("AntiNuMu CC0pi %.1f-%.1f", hAntiNuMuData->GetXaxis()->GetBinLowEdge(x+1), hAntiNuMuData->GetXaxis()->GetBinUpEdge(x+1)));
count++;
}
}
//********************************************************************
void T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint::SetHistograms() {
//********************************************************************
NuMuCC0pi->SetHistograms();
AntiNuMuCC0pi->SetHistograms();
return;
}
//********************************************************************
void T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint::FillHistograms() {
//********************************************************************
NuMuCC0pi->FillHistograms();
AntiNuMuCC0pi->FillHistograms();
return;
}
//********************************************************************
void T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint::ConvertEventRates() {
//********************************************************************
NuMuCC0pi->ConvertEventRates();
AntiNuMuCC0pi->ConvertEventRates();
TH1D* hNuMuCC0pi = (TH1D*)NuMuCC0pi->GetMCHistogram();
TH1D* hAntiNuMuCC0pi = (TH1D*)AntiNuMuCC0pi->GetMCHistogram();
int count = 0;
for (int i = 0; i < hNuMuCC0pi->GetNbinsX(); ++i) {
fMCHist->SetBinContent(count + 1, hNuMuCC0pi->GetBinContent(i + 1));
fMCHist->SetBinError(count + 1, hNuMuCC0pi->GetBinError(i + 1));
count++;
}
for (int i = 0; i < hAntiNuMuCC0pi->GetNbinsX(); ++i) {
fMCHist->SetBinContent(count + 1, hAntiNuMuCC0pi->GetBinContent(i + 1));
fMCHist->SetBinError(count + 1, hAntiNuMuCC0pi->GetBinError(i + 1));
count++;
}
return;
}
diff --git a/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.cxx b/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.cxx
index 0131bc8..bec2fee 100644
--- a/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.cxx
+++ b/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos.cxx
@@ -1,274 +1,274 @@
// 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_NuMu_CC0pi_OC_XSec_2DPcos.h"
static size_t nangbins = 6;
static double angular_binning_costheta[] = {-1, 0., 0.6, 0.75, 0.86, 0.93, 1.};
//********************************************************************
T2K_NuMu_CC0pi_OC_XSec_2DPcos::T2K_NuMu_CC0pi_OC_XSec_2DPcos(nuiskey samplekey) {
//********************************************************************
// Samples overview ---------------------------------------------------
fSettings = LoadSampleSettings(samplekey);
std::string name = fSettings.GetS("name");
std::string descrip = "";
// This has to deal with NuMu on O and on C
if (!name.compare("T2K_NuMu_CC0pi_O_XSec_2DPcos")){
descrip = name +". \n"
"Target: Oxygen \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0pi\n"
- "arXiv:2004.05434";
+ "DOI:10.1103/PhysRevD.101.112004";
fSettings.SetTitle("T2K_NuMu_CC0pi_O_XSec_2DPcos");
fSettings.DefineAllowedTargets("O");
Target = "O";
}
else if (!name.compare("T2K_NuMu_CC0pi_C_XSec_2DPcos")){
descrip = name +". \n"
"Target: Carbon \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0pi\n"
- "arXiv:2004.05434";
+ "DOI:10.1103/PhysRevD.101.112004";
fSettings.SetTitle("T2K_NuMu_CC0pi_C_XSec_2DPcos");
fSettings.DefineAllowedTargets("C");
Target = "C";
}
// Setup common settings
fSettings.SetDescription(descrip);
fSettings.SetXTitle("p_{#mu}-cos#theta_{#mu}");
fSettings.SetYTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)");
fSettings.SetAllowedTypes("DIAG,FULL/FREE,SHAPE,FIX/SYSTCOV/STATCOV","FIX");
fSettings.SetEnuRangeFromFlux(fFluxHist);
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * 1E-38 / (TotalIntegratedFlux()));
// Setup Histograms
SetHistograms();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
bool T2K_NuMu_CC0pi_OC_XSec_2DPcos::isSignal(FitEvent *event){
return SignalDef::isCC0pi(event, 14, EnuMin, EnuMax);
};
void T2K_NuMu_CC0pi_OC_XSec_2DPcos::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_NuMu_CC0pi_OC_XSec_2DPcos::FillHistograms(){
Measurement1D::FillHistograms();
if (Signal){
FillMCSlice( fXVar, fYVar, Weight );
}
}
void T2K_NuMu_CC0pi_OC_XSec_2DPcos::ConvertEventRates(){
for (int i = 0; i < nangbins; i++){
if(Target=="O") fMCHistNuMuO_Slices[i]->GetSumw2();
else if(Target=="C") fMCHistNuMuC_Slices[i]->GetSumw2();
}
// Do standard conversion.
Measurement1D::ConvertEventRates();
// Scale MC slices by their bin width
for (size_t i = 0; i < nangbins; ++i) {
if(Target=="O") fMCHistNuMuO_Slices[i]->Scale(1. / (angular_binning_costheta[i + 1] - angular_binning_costheta[i]));
else if(Target=="C") fMCHistNuMuC_Slices[i]->Scale(1. / (angular_binning_costheta[i + 1] - angular_binning_costheta[i]));
}
// Now Convert into 1D histogram
fMCHist->Reset();
int bincount = 0;
for (int i = 0; i < nangbins; i++){
if(Target=="O"){
for (int j = 0; j < fMCHistNuMuO_Slices[i]->GetNbinsX(); j++){
fMCHist->SetBinContent(bincount+1, fMCHistNuMuO_Slices[i]->GetBinContent(j+1));
bincount++;
}
}
else if(Target=="C"){
for (int j = 0; j < fMCHistNuMuC_Slices[i]->GetNbinsX(); j++){
fMCHist->SetBinContent(bincount+1, fMCHistNuMuC_Slices[i]->GetBinContent(j+1));
bincount++;
}
}
}
return;
}
void T2K_NuMu_CC0pi_OC_XSec_2DPcos::FillMCSlice(double x, double y, double w){
for (size_t i = 0; i < nangbins; ++i) {
if ((y > angular_binning_costheta[i]) && (y <= angular_binning_costheta[i + 1])) {
if(Target=="O") fMCHistNuMuO_Slices[i]->Fill(x, w);
else if(Target=="C") fMCHistNuMuC_Slices[i]->Fill(x, w);
}
}
}
void T2K_NuMu_CC0pi_OC_XSec_2DPcos::SetHistograms(){
// Read covariance matrix
fInputFileCov = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/JointO-C/covmatrix_noreg.root").c_str(),"READ");
TH1D* hLinearResult;
int Nbins;
if(Target=="O"){
fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/JointO-C/linear_unreg_results_O_nuisance.root").c_str(),"READ");
// Read 1D data histogram
hLinearResult = (TH1D*) fInputFile->Get("LinResult");
Nbins = hLinearResult->GetNbinsX();
// Make covariance matrix
fFullCovar = new TMatrixDSym(Nbins);
TMatrixDSym* tempcov = (TMatrixDSym*) fInputFileCov->Get("covmatrixObin");
for(int ibin=0; ibin<Nbins; ibin++) {
for(int jbin=0; jbin<Nbins; jbin++) {
// The factor 1E-2 needed since the covariance matrix in the
// data release is divided by 1E-78
(*fFullCovar)(ibin,jbin) = (*tempcov)(ibin,jbin)*1E-2;
}
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
// Store hLinearResult into fDataHist
fDataHist = new TH1D("LinarResultOxygen","LinarResultOxygen",Nbins,0,Nbins);
for (int bin = 0; bin < Nbins; bin++){
fDataHist->SetBinContent(bin+1, hLinearResult->GetBinContent(bin+1));
}
fDataHist->Reset();
int bincount = 0;
for (int i = 0; i < nangbins; i++){
// Make slices for data
fDataHistNuMuO_Slices.push_back((TH1D*) fInputFile->Get(Form("dataslice_%i",i))->Clone());
fDataHistNuMuO_Slices[i]->SetNameTitle(Form("T2K_NuMu_CC0pi_O_2DPcos_data_Slice%i",i),
(Form("T2K_NuMu_CC0pi_O_2DPcos_data_Slice%i",i)));
fDataHistNuMuO_Slices[i]->Scale(1E-39);
// Loop over nbins and set errors from covariance
for (int j = 0; j < fDataHistNuMuO_Slices[i]->GetNbinsX(); j++){
fDataHistNuMuO_Slices[i]->SetBinError(j+1, sqrt((*fFullCovar)(bincount,bincount))*1E-38);
fDataHist->SetBinContent(bincount+1, fDataHistNuMuO_Slices[i]->GetBinContent(j+1));
fDataHist->SetBinError(bincount+1, fDataHistNuMuO_Slices[i]->GetBinError(j+1));
bincount++;
}
// Save MC slices
fMCHistNuMuO_Slices.push_back((TH1D*) fDataHistNuMuO_Slices[i]->Clone());
fMCHistNuMuO_Slices[i]->SetNameTitle(Form("T2K_NuMu_CC0pi_O_2DPcos_MC_Slice%i",i), (Form("T2K_NuMu_CC0pi_O_2DPcos_MC_Slice%i",i)));
SetAutoProcessTH1(fDataHistNuMuO_Slices[i],kCMD_Write);
SetAutoProcessTH1(fMCHistNuMuO_Slices[i]);
}
}
else if(Target=="C"){
fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/JointO-C/linear_unreg_results_C_nuisance.root").c_str(),"READ");
// Read 1D data histogram
hLinearResult = (TH1D*) fInputFile->Get("LinResult");
Nbins = hLinearResult->GetNbinsX();
// Make the covariance matrix
fFullCovar = new TMatrixDSym(Nbins);
TMatrixDSym* tempcov = (TMatrixDSym*) fInputFileCov->Get("covmatrixCbin");
for(int ibin=0; ibin<Nbins; ibin++) {
for(int jbin=0; jbin<Nbins; jbin++) {
(*fFullCovar)(ibin,jbin) = (*tempcov)(ibin,jbin)*1E-2;
}
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
// Now Convert into 1D histrogram
fDataHist = new TH1D("LinarResultCarbon","LinarResultCarbon",Nbins,0,Nbins);
for (int bin = 0; bin < Nbins; bin++){
fDataHist->SetBinContent(bin+1, hLinearResult->GetBinContent(bin+1));
}
fDataHist->Reset();
int bincount=0;
for (int i = 0; i < nangbins; i++){
// Make slices for data
fDataHistNuMuC_Slices.push_back((TH1D*) fInputFile->Get(Form("dataslice_%i",i))->Clone());
fDataHistNuMuC_Slices[i]->SetNameTitle(Form("T2K_NuMu_CC0pi_C_2DPcos_data_Slice%i",i),
(Form("T2K_NuMu_CC0pi_C_2DPcos_data_Slice%i",i)));
fDataHistNuMuC_Slices[i]->Scale(1E-39);
//Loop over nbins and set errors from covariance
for (int j = 0; j < fDataHistNuMuC_Slices[i]->GetNbinsX(); j++){
fDataHistNuMuC_Slices[i]->SetBinError(j+1, sqrt((*fFullCovar)(bincount,bincount))*1E-38);
fDataHist->SetBinContent(bincount+1, fDataHistNuMuC_Slices[i]->GetBinContent(j+1));
fDataHist->SetBinError(bincount+1, fDataHistNuMuC_Slices[i]->GetBinError(j+1));
bincount++;
}
//Save MC slices
fMCHistNuMuC_Slices.push_back((TH1D*) fDataHistNuMuC_Slices[i]->Clone());
fMCHistNuMuC_Slices[i]->SetNameTitle(Form("T2K_NuMu_CC0pi_C_2DPcos_MC_Slice%i",i), (Form("T2K_NuMu_CC0pi_C_2DPcos_MC_Slice%i",i)));
SetAutoProcessTH1(fDataHistNuMuC_Slices[i],kCMD_Write);
SetAutoProcessTH1(fMCHistNuMuC_Slices[i]);
}
}
return;
};
diff --git a/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.cxx b/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.cxx
index 73644c2..1ad9e6a 100644
--- a/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.cxx
+++ b/src/T2K/T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.cxx
@@ -1,161 +1,161 @@
#include "T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.h"
//********************************************************************
T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint(nuiskey samplekey){
//********************************************************************
fSettings = LoadSampleSettings(samplekey);
std::string descrip = "T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint. \n"
"Target: O-C \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0pi\n"
- "arXiv:2004.05434";
+ "DOI:10.1103/PhysRevD.101.112004";
fSettings.SetTitle("T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("p_{#mu}-cos#theta_{#mu}");
fSettings.SetYTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)");
fSettings.SetAllowedTypes("DIAG,FULL/FREE,SHAPE,FIX/SYSTCOV/STATCOV","FIX");
fSettings.SetEnuRange(0.0, 30.0);
fSettings.DefineAllowedTargets("O,C");
FinaliseSampleSettings();
if (fSubInFiles.size() != 2) {
NUIS_ABORT("T2K NuMuCC0pi O-C joint requires input files in format: NuMuCC0pi on O and NuMuCC0pi on C");
}
std::string inFileNuMuO = fSubInFiles.at(0);
std::string inFileNuMuC = fSubInFiles.at(1);
// Create some config keys
nuiskey NuMuCC0piOKey = Config::CreateKey("sample");
NuMuCC0piOKey.SetS("input", inFileNuMuO);
NuMuCC0piOKey.SetS("type", fSettings.GetS("type"));
NuMuCC0piOKey.SetS("name", "T2K_NuMu_CC0pi_O_XSec_2DPcos");
NuMuCC0piO = new T2K_NuMu_CC0pi_OC_XSec_2DPcos(NuMuCC0piOKey);
nuiskey NuMuCC0piCKey = Config::CreateKey("sample");
NuMuCC0piCKey.SetS("input", inFileNuMuC);
NuMuCC0piCKey.SetS("type", fSettings.GetS("type"));
NuMuCC0piCKey.SetS("name", "T2K_NuMu_CC0pi_C_XSec_2DPcos");
NuMuCC0piC = new T2K_NuMu_CC0pi_OC_XSec_2DPcos(NuMuCC0piCKey);
// Sort out the data hist
this->CombineDataHists();
// Set the covariance
SetCovariance();
// Add to chain for processing
fSubChain.clear();
fSubChain.push_back(NuMuCC0piO);
fSubChain.push_back(NuMuCC0piC);
// This saves information from the sub-measurements
fSaveSubMeas = true;
FinaliseMeasurement();
};
//********************************************************************
void T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::SetCovariance(){
//********************************************************************
// Read covariance matrix
fInputFileCov = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/JointO-C/covmatrix_noreg.root").c_str(),"READ");
TMatrixDSym* tempcov = (TMatrixDSym*) fInputFileCov->Get("covmatrixOeCbin");
fFullCovar = new TMatrixDSym(tempcov->GetNrows());
for(int ibin=0; ibin<tempcov->GetNrows(); ibin++) {
for(int jbin=0; jbin<tempcov->GetNrows(); jbin++) {
// The factor 1E-2 needed since the covariance matrix in the
// data release is divided by 1E-78
(*fFullCovar)(ibin,jbin) = (*tempcov)(ibin,jbin)*1E-2;
}
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
return;
}
//********************************************************************
void T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::CombineDataHists(){
//********************************************************************
TH1D *hNuMuOData = (TH1D*)NuMuCC0piO->GetDataHistogram();
TH1D *hNuMuCData = (TH1D*)NuMuCC0piC->GetDataHistogram();
int nbins = hNuMuOData->GetNbinsX() + hNuMuCData->GetNbinsX();
fDataHist = new TH1D((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str(), nbins, 0, nbins);
fDataHist->SetDirectory(0);
int count = 0;
for (int x=0; x<hNuMuOData->GetNbinsX(); ++x){
fDataHist->SetBinContent(count+1, hNuMuOData->GetBinContent(x+1));
fDataHist->SetBinError(count+1, hNuMuOData->GetBinError(x+1));
fDataHist->GetXaxis()->SetBinLabel(count+1, Form("NuMu CC0pi %.1f-%.1f", hNuMuOData->GetXaxis()->GetBinLowEdge(x+1), hNuMuOData->GetXaxis()->GetBinUpEdge(x+1)));
count++;
}
for (int x=0; x<hNuMuCData->GetNbinsX(); ++x){
fDataHist->SetBinContent(count+1, hNuMuCData->GetBinContent(x+1));
fDataHist->SetBinError(count+1, hNuMuCData->GetBinError(x+1));
fDataHist->GetXaxis()->SetBinLabel(count+1, Form("AntiNuMu CC0pi %.1f-%.1f", hNuMuCData->GetXaxis()->GetBinLowEdge(x+1), hNuMuCData->GetXaxis()->GetBinUpEdge(x+1)));
count++;
}
}
//********************************************************************
void T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::SetHistograms() {
//********************************************************************
NuMuCC0piO->SetHistograms();
NuMuCC0piC->SetHistograms();
return;
}
//********************************************************************
void T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::FillHistograms() {
//********************************************************************
NuMuCC0piO->FillHistograms();
NuMuCC0piC->FillHistograms();
return;
}
//********************************************************************
void T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint::ConvertEventRates() {
//********************************************************************
NuMuCC0piO->ConvertEventRates();
NuMuCC0piC->ConvertEventRates();
TH1D* hNuMuCC0piO = (TH1D*)NuMuCC0piO->GetMCHistogram();
TH1D* hNuMuCC0piC = (TH1D*)NuMuCC0piC->GetMCHistogram();
int count = 0;
for (int i = 0; i < hNuMuCC0piO->GetNbinsX(); ++i) {
fMCHist->SetBinContent(count + 1, hNuMuCC0piO->GetBinContent(i + 1));
fMCHist->SetBinError(count + 1, hNuMuCC0piO->GetBinError(i + 1));
count++;
}
for (int i = 0; i < hNuMuCC0piC->GetNbinsX(); ++i) {
fMCHist->SetBinContent(count + 1, hNuMuCC0piC->GetBinContent(i + 1));
fMCHist->SetBinError(count + 1, hNuMuCC0piC->GetBinError(i + 1));
count++;
}
return;
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sat, May 3, 6:15 AM (1 d, 13 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4971399
Default Alt Text
(32 KB)
Attached To
rNUISANCEGIT nuisancegit
Event Timeline
Log In to Comment