diff --git a/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx b/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx
index 7b0b0d7..e1522bd 100644
--- a/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx
@@ -1,334 +1,226 @@
// 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 "MiniBooNE_NCEL_XSec_Treco_nu.h"
#include "TLorentzVector.h"
-/// ENTIRE CLASS NEEDS FIXING
-
-
// The constructor
MiniBooNE_NCEL_XSec_Treco_nu::MiniBooNE_NCEL_XSec_Treco_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){
-
-
- // // Check if this is a shape only fit - for now, there is no shape option.
- // if (!type.compare("SHAPE")){
- // ERR(WRN) << "MiniBooNE NCEL is not available as a shape only fit... ignoring..." << std::endl;
- // }
+ fName = "MiniBooNE_NCEL_XSec_Treco_nu";
+ fPlotTitles = "; T_{reco} (MeV); Events/(12 MeV)";
+ EnuMin = 0.;
+ EnuMax = 10.0;
+ fTRecoEdges = {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};
- // // Set pointer to the reweighting engine
- // rw_engine = rw;
-
- // // Define the energy of the signal
- // this->EnuMin = 0.;
- // this->EnuMax = 10.;
-
- // // In future read most of these from a card file
- // this->inFile = inputfile;
- // this->fName = "MB_NCEL_XSec_Treco_nu";
- // this->fPlotTitles = "; T_{reco} (MeV); Events/(12 MeV)";
-
- // // Because the binning is in Treco is fairly esoteric, hardcode here
- // this->arr_treco = {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};
-
- // this->SetDataValues(FitPar::GetDataBase()+"/MiniBooNE/ncqe/input_data.txt");
- // this->SetCovarMatrix(FitPar::GetDataBase()+"/MiniBooNE/ncqe/ErrorMatrix.tab", 51);
- // this->SetResponseMatrix(FitPar::GetDataBase()+"/MiniBooNE/ncqe/response_mat.txt", 51);
-
- // // Check if we're using fake data
- // if (!fakeDataFile.empty()) this->SetFakeDataValues(fakeDataFile);
-
- // // This will be the final data histogram
- // this->fMCHist = new TH1D((this->fName+"_MC").c_str(), (this->fName+this->fPlotTitles).c_str(), 51, this->arr_treco);
-
- // // Usually, the fMCFine histogram is a finer binned version of fMCHist. But as NCEL requires a Ttrue histogram, co-opt fMCFine for this purpose.
- // // Should probably change the naming to reflect the possible other use of this histogram.
- // this->fMCFine = new TH1D((this->fName+"_Ttrue").c_str(), (this->fName+this->fPlotTitles).c_str(), 50, 0, 900);
-
- // // Read in the histograms from the NEUT file that are required for normalisation
- // TFile *in = new TFile(this->inFile.c_str());
- // this->fFluxHist = (TH1D*)in->Get((PlotUtils::GetObjectWithName(in, "flux")).c_str());
- // this->fFluxHist->SetNameTitle((this->fName+"_FLUX").c_str(), (this->fName+";E_{#nu} (GeV)").c_str());
-
- // this->fEventHist = (TH1D*)in->Get((PlotUtils::GetObjectWithName(in, "evtrt")).c_str());
- // this->fEventHist->SetNameTitle((this->fName+"_EVT").c_str(), (this->fName+";E_{#nu} (GeV); Event Rate").c_str());
-
- // // Read in the file once only
- // tn = new TChain("neuttree", "");
- // tn->Add(Form("%s/neuttree", this->inFile.c_str()));
- // fNEvents = tn->GetEntries();
- // nvect = NULL;
- // tn->SetBranchAddress("vectorbranch", &nvect);
-
- // // 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.
- // double nMolMB = 6.023E+23*0.845*4.0*M_PI*610.6*610.6*610.6/3.0;
- // this->fScaleFactor = (this->fEventHist->Integral()*1E-38*14.08/(fNEvents+0.))*nMolMB*0.646165;
-};
+ SetDataValues(FitPar::GetDataBase()+"/MiniBooNE/ncqe/input_data.txt");
+ SetCovarMatrix(FitPar::GetDataBase()+"/MiniBooNE/ncqe/ErrorMatrix.tab", 51);
+ SetResponseMatrix(FitPar::GetDataBase()+"/MiniBooNE/ncqe/response_mat.txt", 51);
+ // Setup MC Hists
+ Measurement1D::SetupDefaultHist();
-void MiniBooNE_NCEL_XSec_Treco_nu::Reconfigure(double norm, bool fullconfig){
-
- // // Clear the current histogram before repopulating
- // this->fMCHist->Reset();
- // this->fMCFine->Reset();
- // this->fCurrentNorm = norm;
-
- // // Loop over all events at each iteration of the fit
- // for (int i = 0; i < fNEvents; ++i){
- // tn->GetEntry(i);
-
- // if (!isSignal(nvect)) continue;
-
- // // Find the weight
- // double rw_weight = 1;
- // rw_weight = rw_engine->CalcWeight(customEvent);
-
- // // Skip any events with suspiciously large weights...
- // if (rw_weight > 200) {ERR(WRN) << "LARGE WEIGHT: " << rw_weight << std::endl; break;}
-
- // // Sum of the true kinetic energies of particles
- // double t_true = 0.;
-
- // // Loop over the particle stack
- // for (UInt_t j = 2; j < nvect->Npart(); ++j){
-
- // // Add the kinetic energies of any nucleons
- // if (abs((nvect->PartInfo(j))->fPID) == 2212)
- // t_true += (nvect->PartInfo(j))->fP.E()/1000 - PhysConst::mass_proton;
- // else if (abs((nvect->PartInfo(j))->fPID) == 2112)
- // t_true += (nvect->PartInfo(j))->fP.E()/1000 - PhysConst::mass_neutron;
-
- // }
-
- // // Now fill the Ttrue histogram
- // this->fMCFine->Fill(t_true*1000., rw_weight);
- // }
+ // 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);
- // // 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)*this->response_mat->GetBinContent(ttrue+1, treco+1);
- // this->fMCHist->SetBinContent(treco+1, total);
- // }
+ // 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.
+ double nMolMB = 6.023E+23*0.845*4.0*M_PI*610.6*610.6*610.6/3.0;
+ fScaleFactor = (fEventHist->Integral("width")*1E-38*14.08/(fNEvents+0.))*nMolMB*0.646165;
- // // Scale
- // this->fMCHist->Scale(this->fScaleFactor, "width");
- // this->fMCFine->Scale(this->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);
- // }
+ void MiniBooNE_NCEL_XSec_Treco_nu::FillEventVariables(FitEvent *event){
- // // Normalisation factor if one has been provided.
- // if (norm){
- // this->fMCHist->Scale(1.0/norm);
- // } else {
- // this->fMCHist->Scale(0);
- // }
-
- // return;
-};
+ 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;
+}
-bool MiniBooNE_NCEL_XSec_Treco_nu::isSignal(NeutVect *nvect){
-
- // Only interested in true NCEL events
- // if (nvect->Mode != 51 && nvect->Mode != 52) return false;
+bool MiniBooNE_NCEL_XSec_Treco_nu::ScaleEvents(){
- // // Only look at numu events
- // if ((nvect->PartInfo(0))->fPID != 14 && (nvect->PartInfo(0))->fPID != 12) return false;
+ // 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);
+ }
- // // Restrict energy range
- // if ((nvect->PartInfo(0))->fP.E() < this->EnuMin*1000 || (nvect->PartInfo(0))->fP.E() > this->EnuMax*1000) return false;
+ // 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;
};
// 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);
-
+
// 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;
};
-// The covariance matrix contains all of the information for the chi2 calculation
-double MiniBooNE_NCEL_XSec_Treco_nu::GetChi2(){
- // double chi2 = 0;
-
- // int nBins = this->fDataHist->GetNbinsX();
-
- // for (int i = 0; i < nBins; ++i){
- // for (UInt_t j = 0; j < nBins; ++j){
-
- // double iDiff = this->fDataHist->GetBinContent(i+1) - this->fMCHist->GetBinContent(i+1);
- // double jDiff = this->fDataHist->GetBinContent(j+1) - this->fMCHist->GetBinContent(j+1);
- // chi2 += iDiff*(*this->covar)(i, j)*jDiff;
- // }
- // }
- // return chi2;
-
-};
-
-
-// Function to make an Asimov dataset (so don't throw any errors)
-void MiniBooNE_NCEL_XSec_Treco_nu::SetFakeDataValues(std::string fakeDataFile) {
-
- // // This is the published data
- // TH1D *tempData = (TH1D*)this->fDataHist->Clone();
- // TFile *fake = new TFile(fakeDataFile.c_str());
-
- // // This is the fake data
- // this->fDataHist = (TH1D*)fake->Get((this->fName+"_MC").c_str());
- // this->fDataHist ->SetNameTitle((this->fName+"_FAKE").c_str(), (this->fName+this->fPlotTitles).c_str());
-
- // for (int xBin = 0; xBin < this->fDataHist->GetNbinsX(); ++xBin){
-
- // // If the fake data or real didn't didn't fill the bin, can't assign an error
- // if (!this->fDataHist->GetBinContent(xBin+1) || !tempData->GetBinContent(xBin+1)){
- // this->fDataHist->SetBinError(xBin+1, 0);
- // continue;
- // }
-
- // double err = tempData->GetBinError(xBin+1)*
- // this->fDataHist->GetBinContent(xBin+1)/(tempData->GetBinContent(xBin+1)+0.);
- // this->fDataHist->SetBinError(xBin+1, err);
- // }
-
- // delete tempData;
- // 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){
// std::string line;
// std::ifstream input(inputFile.c_str(),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, this->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){
// Make a counter to track the line number
// int xBin = 0;
// std::string line;
// std::ifstream response(responseFile.c_str(),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, this->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_NCEL_XSec_Treco_nu.h b/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.h
index 33b3b09..5885153 100644
--- a/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.h
+++ b/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.h
@@ -1,71 +1,56 @@
// 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/>.
*******************************************************************************/
#ifndef MINIBOONE_NCEL_XSEC_TRECO_NU_H_SEEN
#define MINIBOONE_NCEL_XSEC_RECO_NU_H_SEEN
#include "Measurement1D.h"
//********************************************************************
class MiniBooNE_NCEL_XSec_Treco_nu : public Measurement1D {
//********************************************************************
public:
MiniBooNE_NCEL_XSec_Treco_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile);
virtual ~MiniBooNE_NCEL_XSec_Treco_nu() {};
- void Reconfigure(double norm, bool fullconfig=false);
+ void FillEventVariables(FitEvent *event);
- double GetChi2();
+ bool isSignal(FitEvent *event);
-#ifdef __NEUT_ENABLED__
- bool isSignal(NeutVect *nvect);
-#endif
+ void ScaleEvents();
void SetDataValues(std::string inputFile);
- void SetFakeDataValues(std::string fakeDataFile);
-
private:
void SetCovarMatrix(std::string covarFile, int dim);
void SetResponseMatrix(std::string responseFile, int dim);
- std::string inFile;
- TChain *tn;
- Int_t fNEvents;
-#ifdef __NEUT_ENABLED__
- NeutVect *nvect;
-#endif
-
// Because the Treco binning is irregular, store an array of bin edges...
double arr_treco[52];
TH1D *BKGD_other;
TH1D *BKGD_irrid;
TH2D *response_mat;
- Double_t totIntFlux;
-
- FitWeight *rw_engine;
-
};
#endif