diff --git a/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx b/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx
index e46ae68..fc704e7 100644
--- a/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx
@@ -1,264 +1,264 @@
// 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 "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};
-
+ 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();
// 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");
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");
this->newFluxHist->Print();
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(),ifstream::in);
+ std::ifstream covar(covarFile.c_str(),std::ifstream::in);
this->covar = new TMatrixDSym(dim);
if(covar.is_open()) LOG(DEB) << "Reading covariance matrix from file: " << covarFile << std::endl;
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(),ifstream::in);
+ std::ifstream input(inputFile.c_str(),std::ifstream::in);
if(input.is_open()) LOG(DEB) << "Reading data from file: " << inputFile << std::endl;
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(),ifstream::in);
+ 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()) LOG(DEB) << "Reading in the response matrix from file: " << responseFile << std::endl;
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 4c24717..a69e496 100644
--- a/src/MiniBooNE/MiniBooNE_NCpi0_XSec_1Dppi0_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_NCpi0_XSec_1Dppi0_nu.cxx
@@ -1,171 +1,170 @@
// 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 "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) {
LOG(SAM) << this->fName << "Setting data for " << this->fName << std::endl;
LOG(SAM) << this->fName << "From: " << dataFile << std::endl;
LOG(SAM) << this->fName << "Reading error from covariance" << std::endl;
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) {
LOG(SAM) << this->fName << "===============" << std::endl;
LOG(SAM) << this->fName << "Reading covariance: " << this->fName << std::endl;
LOG(SAM) << this->fName << "From: " << covarFile << std::endl;
// tracks line number
int row = 0;
std::string line;
- std::ifstream covar(covarFile.c_str(), ifstream::in);
+ std::ifstream covar(covarFile.c_str(), std::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;
};
-