diff --git a/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTpi_nu.cxx b/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTpi_nu.cxx
index 369573c..e591cf4 100644
--- a/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTpi_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTpi_nu.cxx
@@ -1,90 +1,90 @@
// 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 "MiniBooNE_CC1pip_XSec_1DTpi_nu.h"
//********************************************************************
MiniBooNE_CC1pip_XSec_1DTpi_nu::MiniBooNE_CC1pip_XSec_1DTpi_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_CC1pip_XSec_1DTpi_nu sample. \n" \
"Target: CH \n" \
"Flux: MiniBooNE Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 muon, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("T_{#pi} (MeV)");
- fSettings.SetYTitle("d#sigma/dT_{#pi^{+}} (cm^{2}/MeV/CH_{2})");
+ fSettings.SetYTitle("d#sigma/dT_{#pi} (cm^{2}/MeV/CH_{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.SetNormError(0.107);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MiniBooNE_CC1pip_XSec_1DTpi_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "MiniBooNE/CC1pip/ccpipXSec_KEpi.txt" );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width")*double(1E-38)/double(fNEvents)*(14.08)/TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromDiagonal();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MiniBooNE_CC1pip_XSec_1DTpi_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
if (event->NumFSParticle(211) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
// No W cut on MiniBooNE data from publication
// WARNING: DRAGONS LAY HERE! Mike Wilking's thesis might not have this. Beware that the publication says W < 1.35 GeV, but this is "efficiency corrected", i.e. FILLED WITH MONTE-CARLO!!!!!!!! AAAAAH
//double hadMass = FitUtils::Wrec(Pnu, Pmu, Ppip);
double Tpi = FitUtils::T(Ppip)*1000.;
fXVar = Tpi;
return;
};
bool MiniBooNE_CC1pip_XSec_1DTpi_nu::isSignal(FitEvent *event) {
return SignalDef::isCC1pi(event, 14, 211, EnuMin, EnuMax);
}
diff --git a/src/MiniBooNE/MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.cxx b/src/MiniBooNE/MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.cxx
index 50d934e..a4e0f51 100644
--- a/src/MiniBooNE/MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.cxx
@@ -1,85 +1,85 @@
// 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 "MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.h"
//********************************************************************
MiniBooNE_CC1pip_XSec_2DQ2Enu_nu::MiniBooNE_CC1pip_XSec_2DQ2Enu_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_CC1pip_XSec_2DQ2Enu_nu sample. \n" \
"Target: CH \n" \
"Flux: MiniBooNE Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
- fSettings.SetXTitle("E_{#nu} (MeV); Q^{2} (MeV^{2}/c^{4})");
- fSettings.SetYTitle("Q^{2} (MeV^{2}/c^{4})");
- fSettings.SetZTitle(" d#sigma(E_{#nu})/dQ^{2} (cm^{2}/(MeV^{2}/c^{4})/CH_{2})");
+ fSettings.SetXTitle("E_{#nu} (MeV)");
+ fSettings.SetYTitle("Q^{2} (MeV^{2})");
+ fSettings.SetZTitle("d#sigma(E_{#nu})/dQ^{2} (cm^{2}/MeV^{2}/CH_{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 3.0);
fSettings.SetNormError(0.107);
fSettings.DefineAllowedTargets("C,H");
fSettings.SetTitle("MiniBooNE_CC1pip_XSec_2DQ2Enu_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MiniBooNE/CC1pip/ccpipXSecs.root;QSQVENUXSec" );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents) * (14.08);
// Plot Setup -------------------------------------------------------
SetDataFromRootFile( fSettings.GetDataInput() );
SetCovarFromDiagonal();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MiniBooNE_CC1pip_XSec_2DQ2Enu_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
if (event->NumFSParticle(13) == 0) return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double Enu = Pnu.E();
double Q2 = -1*(Pnu-Pmu).Mag2();
fXVar = Enu;
fYVar = Q2;
return;
};
//********************************************************************
bool MiniBooNE_CC1pip_XSec_2DQ2Enu_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCC1pi(event, 14, 211, EnuMin, EnuMax);
}
diff --git a/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DEnu_nu.cxx b/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DEnu_nu.cxx
index 16e6a56..1624355 100644
--- a/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DEnu_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DEnu_nu.cxx
@@ -1,87 +1,87 @@
// 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 "MiniBooNE_CCQE_XSec_1DEnu_nu.h"
//********************************************************************
MiniBooNE_CCQE_XSec_1DEnu_nu::MiniBooNE_CCQE_XSec_1DEnu_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_CCQE_XSec_1DEnu_nu sample. \n"
"Target: CH2 \n"
"Flux: MiniBooNE Forward Horn Current\n"
"Signal: Any event with 1 muon, any nucleons, and no "
"other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("E_{#nu} (GeV)");
- fSettings.SetYTitle("#sigma(E_{#nu}) (cm^{2}/Neutron)");
+ fSettings.SetYTitle("#sigma(E_{#nu}) (cm^{2}/neutron)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.4, 2.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.FoundFill("name", "CCQELike", ccqelike, true);
// CCQELike plot information
fSettings.SetTitle("MiniBooNE_CCQE_XSec_1DEnu_nu");
if (ccqelike) {
fSettings.SetDataInput(FitPar::GetDataBase() +
"MiniBooNE/ccqe/asne_like.txt");
} else {
fSettings.SetDataInput(FitPar::GetDataBase() +
"MiniBooNE/ccqe/asne_con.txt");
}
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2
fScaleFactor =
GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents) * (14.08/6.0);
// Plot Setup -------------------------------------------------------
SetDataFromTextFile(fSettings.GetDataInput());
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void MiniBooNE_CCQE_XSec_1DEnu_nu::FillEventVariables(FitEvent *event) {
if (isSignal(event)) {
fXVar = event->GetNeutrinoIn()->fP.E() * 1E-3;
}
};
bool MiniBooNE_CCQE_XSec_1DEnu_nu::isSignal(FitEvent *event) {
// If CC0pi, include both charges
if (ccqelike) {
if (SignalDef::isCC0pi(event, 14, EnuMin, EnuMax) ||
SignalDef::isCC0pi(event, -14, EnuMin, EnuMax))
return true;
} else {
if (SignalDef::isCCQELike(event, 14, EnuMin, EnuMax))
return true;
}
return false;
}
diff --git a/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx b/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx
index 33390da..be6f4cf 100644
--- a/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx
@@ -1,293 +1,295 @@
// 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 "MiniBooNE_NCEL_XSec_Treco_nu.h"
#include "TLorentzVector.h"
//********************************************************************
MiniBooNE_NCEL_XSec_Treco_nu::MiniBooNE_NCEL_XSec_Treco_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_NCEL_XSec_Treco_nu sample. \n"
"Target: CH \n"
"Flux: MiniBooNE Numu Flux \n"
"Signal: Any event with True NCEL modes \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("T_{reco} (MeV)");
fSettings.SetYTitle("Events/(12 MeV)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/NORM/MASK",
"FIX/FULL,DIAG");
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetEnuRange(0.0, 10.0);
double arr_treco[52] = {
40.0, 52.0, 63.9, 75.9, 87.8, 99.8, 111.8, 123.7, 135.7,
147.6, 159.6, 171.6, 183.5, 195.5, 207.5, 219.4, 231.4, 243.3,
255.3, 267.3, 279.2, 291.2, 303.1, 315.1, 327.1, 339.0, 351.0,
362.9, 374.9, 386.9, 398.8, 410.8, 422.7, 434.7, 446.7, 458.6,
470.6, 482.5, 494.5, 506.5, 518.4, 530.4, 542.4, 554.3, 566.3,
578.2, 590.2, 602.2, 614.1, 626.1, 638.0, 650.0};
SetDataValues(FitPar::GetDataBase() + "/MiniBooNE/ncqe/input_data.txt",
arr_treco);
SetCovarMatrix(FitPar::GetDataBase() + "/MiniBooNE/ncqe/ErrorMatrix.tab", 51);
SetResponseMatrix(FitPar::GetDataBase() + "/MiniBooNE/ncqe/response_mat.txt",
51, arr_treco);
SetFluxHistogram(FitPar::GetDataBase() + "/MiniBooNE/ncqe/flux.txt");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// The scale factor is quite complicated because MB didn't divide by number of
// targets. nMolMB is the number of CH_2 molecules in the MB FV (610.6 cm
// radius sphere) and 0.845 is the published density of the mineral oil.
// UPDATE: nMolMB is the number of CH_2 molecules in the MB FV (500 cm radius
// sphere) and 0.845 is the published density of the mineral oil. UPDATE
// UPDATE: They didn't account for fiducial cut, so neither do we.
double nMolMB = 6.023E+23 * 0.845 * 4.0 * M_PI * 610.6 * 610.6 * 610.6 / 3.0;
double POT = 6.46e20;
// Need to update scalefactor to reflect actual flux used
fScaleFactor = (this->newFluxHist->Integral("") * POT *
(GetEventHistogram()->Integral("width")) * 1E-38 * 14.08 /
(fNEvents + 0.)) *
nMolMB;
fScaleFactor /= GetFluxHistogram()->Integral("width");
// Final setup ---------------------------------------------------
FinaliseMeasurement();
+ this->fPlotTitles = fSettings.GetFullTitles();
+
// Usually the MCFine histogram is a finer binned version of MC Hist.
// In this case we need to use it to save the true distribution before
// smearing.
if (fMCFine)
delete fMCFine;
fMCFine = new TH1D((this->fName + "_Ttrue").c_str(),
(this->fName + this->fPlotTitles).c_str(), 50, 0, 900);
};
void MiniBooNE_NCEL_XSec_Treco_nu::Write(std::string arg) {
- newFluxHist->Write("MB_NCEL_newFlux");
- response_mat->Write("MB_NCEL_response matrix");
+ // newFluxHist->Write("MB_NCEL_newFlux");
+ response_mat->Write("MB_NCEL_response_matrix");
Measurement1D::Write(arg);
return;
}
void MiniBooNE_NCEL_XSec_Treco_nu::FillEventVariables(FitEvent *event) {
double t_raw = 0.0;
// Loop and add all Tnucleon
for (UInt_t i = 0; i < event->Npart(); i++) {
if (event->PartInfo(i)->Status() != kFinalState)
continue;
int pdg = event->PartInfo(i)->fPID;
if (pdg == 2212 || pdg == 2112) {
t_raw += FitUtils::T(event->PartInfo(i)->fP) * 1.E3;
}
}
fXVar = t_raw;
}
void MiniBooNE_NCEL_XSec_Treco_nu::ScaleEvents() {
// Now convert Ttrue to Treco...
for (int treco = 0; treco < 51; ++treco) {
double total = 0.;
for (int ttrue = 0; ttrue < 50; ++ttrue)
total += fMCFine->GetBinContent(ttrue + 1) *
response_mat->GetBinContent(ttrue + 1, treco + 1);
fMCHist->SetBinContent(treco + 1, total);
}
// Scale
this->fMCHist->Scale(this->fScaleFactor, "width");
this->fMCFine->Scale(this->fScaleFactor, "width");
PlotUtils::ScaleNeutModeArray((TH1 **)fMCHist_PDG, fScaleFactor, "width");
// Add in the backgrounds...
for (int treco = 0; treco < 51; ++treco) {
double total = this->fMCHist->GetBinContent(treco + 1) +
this->BKGD_other->GetBinContent(treco + 1) +
this->BKGD_irrid->GetBinContent(treco + 1);
this->fMCHist->SetBinContent(treco + 1, total);
}
}
bool MiniBooNE_NCEL_XSec_Treco_nu::isSignal(FitEvent *event) {
// Should put in MB SignalDef eventually
if (event->Mode != 51 && event->Mode != 52)
return false;
// Numu or nue
if (event->PDGnu() != 14 && event->PDGnu() != 12)
return false;
// Enu
if (event->Enu() < EnuMin * 1000.0 || event->Enu() > EnuMax * 1000.0)
return false;
return true;
};
void MiniBooNE_NCEL_XSec_Treco_nu::SetFluxHistogram(std::string dataFile) {
this->newFluxHist = PlotUtils::GetTH1DFromFile(
dataFile.c_str(), (this->fName + "Real Flux Hist"), "idgaf");
return;
}
// Read in the covariance matrix from the file specified in the constructor
void MiniBooNE_NCEL_XSec_Treco_nu::SetCovarMatrix(std::string covarFile,
int dim) {
// Use Utils
// // Make a counter to track the line number
int row = 0;
std::string line;
std::ifstream covar(covarFile.c_str(), std::ifstream::in);
this->covar = new TMatrixDSym(dim);
if (covar.is_open())
NUIS_LOG(DEB, "Reading covariance matrix from file: " << covarFile);
while (std::getline(covar >> std::ws, line, '\n')) {
std::istringstream stream(line);
double entry;
int column = 0;
// Loop over entries and insert them into matrix
// Multiply by the errors to get the covariance, rather than the correlation
// matrix
while (stream >> entry) {
(*this->covar)(row, column) = entry;
if (row == column)
this->fDataHist->SetBinError(row + 1, sqrt(entry));
column++;
}
row++;
}
// // Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
// Override the usual function in the base class because this is more
// complicated for the NCEL sample...
void MiniBooNE_NCEL_XSec_Treco_nu::SetDataValues(std::string inputFile,
double *arr_treco) {
std::string line;
std::ifstream input(inputFile.c_str(), std::ifstream::in);
if (input.is_open()) {
NUIS_LOG(DEB, "Reading data from file: " << inputFile);
}
this->fDataHist =
new TH1D((this->fName + "_data").c_str(),
(this->fName + this->fPlotTitles).c_str(), 51, arr_treco);
this->BKGD_other =
new TH1D((this->fName + "_BKGD_other").c_str(),
(this->fName + this->fPlotTitles).c_str(), 51, arr_treco);
this->BKGD_irrid =
new TH1D((this->fName + "_BKGD_irrid").c_str(),
(this->fName + this->fPlotTitles).c_str(), 51, arr_treco);
// To get the nDOF correct...
this->fNDataPointsX = 52;
double entry = 0;
int xBin = 0;
// First line is the MB data
std::getline(input >> std::ws, line, '\n');
std::istringstream stream1(line);
while (stream1 >> entry) {
this->fDataHist->SetBinContent(xBin + 1, entry);
xBin++;
}
// Second line is "other" backgrounds
std::getline(input >> std::ws, line, '\n');
std::istringstream stream2(line);
entry = 0;
xBin = 0;
while (stream2 >> entry) {
this->BKGD_other->SetBinContent(xBin + 1, entry);
xBin++;
}
// Third line is the irreducible background
std::getline(input >> std::ws, line, '\n');
std::istringstream stream3(line);
entry = 0;
xBin = 0;
while (stream3 >> entry) {
this->BKGD_irrid->SetBinContent(xBin + 1, entry);
xBin++;
}
};
// Read in the response matrix -- thus far, a response matrix is unique to the
// NCEL sample
void MiniBooNE_NCEL_XSec_Treco_nu::SetResponseMatrix(std::string responseFile,
int dim,
double *arr_treco) {
// Make a counter to track the line number
int xBin = 0;
std::string line;
std::ifstream response(responseFile.c_str(), std::ifstream::in);
// Response matrix: x axis is Ttrue, y axis is Treco
this->response_mat = new TH2D((this->fName + "_RESPONSE_MATRIX").c_str(),
(this->fName + this->fPlotTitles).c_str(), 50,
0, 900, 51, arr_treco);
if (response.is_open()) {
NUIS_LOG(DEB, "Reading in the response matrix from file: " << responseFile);
}
while (std::getline(response, line, '\n')) {
std::istringstream stream(line);
double entry;
int yBin = 0;
// Loop over entries and insert them into matrix
// Multiply by the errors to get the covariance, rather than the correlation
// matrix
while (stream >> entry) {
this->response_mat->SetBinContent(xBin + 1, yBin + 1, entry);
yBin++;
}
xBin++;
}
};
diff --git a/src/MiniBooNE/MiniBooNE_NCpi0_XSec_1Dppi0_nu.cxx b/src/MiniBooNE/MiniBooNE_NCpi0_XSec_1Dppi0_nu.cxx
index 45cd717..809d68a 100644
--- a/src/MiniBooNE/MiniBooNE_NCpi0_XSec_1Dppi0_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_NCpi0_XSec_1Dppi0_nu.cxx
@@ -1,188 +1,146 @@
// 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 "MiniBooNE_NCpi0_XSec_1Dppi0_nu.h"
// The constructor
MiniBooNE_NCpi0_XSec_1Dppi0_nu::MiniBooNE_NCpi0_XSec_1Dppi0_nu(
std::string inputfile, FitWeight *rw, std::string type,
std::string fakeDataFile, Double_t *anuBins)
: isComb(false) {
if (anuBins != NULL)
isComb = true;
isComb = false;
- // Needs Updating
-
- // // Set pointer to the reweighting engine
- // rw_engine = rw;
- // this->fBeamDistance = 0.541;
-
- // // Define the energy region
- // this->EnuMin = 0.;
- // this->EnuMax = 4.;
-
- // // In future read most of these from a card file
- // this->inFile = inputfile;
- // this->fName = "MB_NCpi0_XSec_numu_1Dppi0";
- // this->fPlotTitles = "; p_{#pi^{0}} (GeV/c); d#sigma/dp_{#pi^{0}}
- // (cm^{2}/(GeV/c)/nucleon)";
- // this->SetCovarMatrix(FitPar::GetDataBase()+"/MiniBooNE/nc1pi0/nuppi0xsecerrormatrix.txt",
- // 11);
- // this->SetDataValues(FitPar::GetDataBase()+"/MiniBooNE/nc1pi0/nuppi0xsec_edit.txt");
- // this->fNormError=0.107;
-
- // if (isComb) {
- // fName += "_comb";
- // this->fNDataPointsX = 11;
- // this->fXBins = anuBins;
- // }
-
- // this->fMCHist = new TH1D((this->fName+"_MC").c_str(),
- // (this->fName+this->fPlotTitles).c_str(), this->fNDataPointsX-1,
- // this->fXBins); this->fMCFine = new TH1D((this->fName+"_MC_FINE").c_str(),
- // (this->fName+this->fPlotTitles).c_str(), (this->fNDataPointsX - 1)*10,
- // this->fXBins[0], this->fXBins[this->fNDataPointsX -1]);
-
- // this->ReadEventFile();
-
- // // Different generators require slightly different rescaling factors.
- // if (this->fEventType == 0) this->fScaleFactor =
- // (GetEventHistogram()->Integral("width")*1E-38/(fNEvents+0.))*14.08/14.0/this->TotalIntegratedFlux();
- // // NEUT else if (this->fEventType == 1) this->fScaleFactor =
- // (GetEventHistogram()->Integral()*1E-38/(fNEvents+0.))*14.08*6.0/14./GetFluxHistogram()->Integral();
- // // NUWRO else if (this->fEventType == 5) this->fScaleFactor =
- // (GetEventHistogram()->Integral()*1E-38/(fNEvents+0.))*14.08*6.0/14./GetFluxHistogram()->Integral();
- // // GENIE
};
void MiniBooNE_NCpi0_XSec_1Dppi0_nu::FillEventVariables(FitEvent *event) {
TLorentzVector Pnu = (event->PartInfo(0))->fP;
TLorentzVector Pmu;
TLorentzVector Ppi0;
double EHad = 0;
pi0Cnt = 0;
bad_particle = false;
for (UInt_t j = 2; j < event->Npart(); ++j) {
if (!((event->PartInfo(j))->fIsAlive) &&
(event->PartInfo(j))->fNEUTStatusCode != 0)
continue;
int PID = (event->PartInfo(j))->fPID;
double KE = (event->PartInfo(j))->fP.E() - (event->PartInfo(j))->fMass;
if (PID == 111) {
Ppi0 = event->PartInfo(j)->fP;
EHad += KE;
} else if (PID == 2112 || PID == 2212)
EHad += KE;
else if (PID == -13)
Pmu = event->PartInfo(j)->fP;
if (abs(PID) >= 113 && abs(PID) <= 557)
bad_particle = true;
else if (abs(PID) == 11 || abs(PID) == 13 || abs(PID) == 15 ||
abs(PID) == 17)
bad_particle = true;
else if (PID == 111)
pi0Cnt++;
}
double bind = 34.0;
if (isComb)
bind = 30.0;
// double hadMass = FitUtils::Wrec(Pnu, Pmu, Ppi0);
double ppi0 = Ppi0.Vect().Mag() / 1000.0;
fXVar = ppi0;
return;
};
bool MiniBooNE_NCpi0_XSec_1Dppi0_nu::isSignal(FitEvent *event) {
return SignalDef::isNC1pi(event, 14, 111, EnuMin, EnuMax);
};
void MiniBooNE_NCpi0_XSec_1Dppi0_nu::SetDataValues(std::string dataFile) {
NUIS_LOG(SAM, this->fName << "Setting data for " << this->fName);
NUIS_LOG(SAM, this->fName << "From: " << dataFile);
NUIS_LOG(SAM, this->fName << "Reading error from covariance");
TGraph *gr = new TGraph(dataFile.c_str());
this->fXBins = gr->GetX();
this->fDataValues = gr->GetY();
this->fNDataPointsX = gr->GetN();
// get the diagonal elements
int rows = (this->tempCovar)->GetNrows();
Double_t errors[rows + 1];
for (int i = 0; i < rows; i++)
errors[i] = sqrt((*this->tempCovar)(i, i) * 1E-81);
errors[rows] = 0.;
this->fDataErrors = errors;
this->fDataHist = new TH1D((this->fName + "_data").c_str(),
(this->fName + this->fPlotTitles).c_str(),
this->fNDataPointsX - 1, this->fXBins);
for (int i = 0; i < this->fNDataPointsX; ++i) {
this->fDataHist->SetBinContent(i + 1, this->fDataValues[i]);
this->fDataHist->SetBinError(i + 1, this->fDataErrors[i]);
}
return;
}
void MiniBooNE_NCpi0_XSec_1Dppi0_nu::SetCovarMatrix(std::string covarFile,
int dim) {
NUIS_LOG(SAM, this->fName << "===============");
NUIS_LOG(SAM, this->fName << "Reading covariance: " << this->fName);
NUIS_LOG(SAM, this->fName << "From: " << covarFile);
// tracks line number
int row = 0;
std::string line;
std::ifstream covar(covarFile.c_str(), ifstream::in);
this->tempCovar = new TMatrixDSym(dim);
// while we're on a line in covar
while (std::getline(covar >> std::ws, line, '\n')) {
std::istringstream stream(line);
// this is the netry we're reading!
double entry;
// this is the column counter!
int column = 0;
while (stream >> entry) {
// get the covariance entry.
// 1E-81 from the data release listing this unit
double val = entry; // * 1E-81;
// then fill the covariance matrix's row and column with this value,
(*this->tempCovar)(row, column) = val;
column++;
}
row++;
}
this->covar = (TMatrixDSym *)this->tempCovar->Clone();
TDecompChol LU = TDecompChol(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
(*this->covar) *= 1E81 * 1E-76;
return;
};