diff --git a/data/MINERvA/CC0pi_2D/minerva_cc0pi_ME_ptpz.root b/data/MINERvA/CC0pi_2D/minerva_cc0pi_ME_ptpz.root
new file mode 100644
index 0000000..0d05216
Binary files /dev/null and b/data/MINERvA/CC0pi_2D/minerva_cc0pi_ME_ptpz.root differ
diff --git a/data/MINERvA/CC0pi_2D_antinu/data_release_pzpt.root b/data/MINERvA/CC0pi_2D_antinu/data_release_pzpt.root
new file mode 100644
index 0000000..f6f640a
Binary files /dev/null and b/data/MINERvA/CC0pi_2D_antinu/data_release_pzpt.root differ
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2D_antinu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_2D_antinu.cxx
index 4074216..bb5980d 100644
--- a/src/MINERvA/MINERvA_CC0pi_XSec_2D_antinu.cxx
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_2D_antinu.cxx
@@ -1,199 +1,226 @@
-// Adrian Orea
-// I used the file MINERvA_CCinc_XSec_2DEavq3_nu.cxx as a template
-// Also, I am fully aware of the naming typo (should be ptpz), but Everything is already named the same way so...
-
// 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 .
*******************************************************************************/
/*
Author : Adrian Orea
Clarence Wret Cleanup 2019: Data was missing
Signal definition wrong (missing 120 MeV KE cut, missing pz and pt cut)
Not fully implemented
Assumed generator send neutrinos along z-axis
*/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CC0pi_XSec_2D_antinu.h"
//********************************************************************
void MINERvA_CC0pi_XSec_2D_antinu::SetupDataSettings() {
//********************************************************************
+ if (IsME) {
+ std::cerr << "Medium energy implementation does not support covariance" << std::endl;
+ std::cerr << "If you want to run anyway, comment me out" << std::endl;
+ std::cerr << "I'm in " << __FILE__ << ":" << __LINE__ << std::endl;
+ throw;
+ }
+
// Set Distribution
// See header file for enum and some descriptions
std::string name = fSettings.GetS("name");
- // We're lucky to have three different MINERvA CC0pi anti-numu 2D distributions
- if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_antinu")) fDist = kPtPz;
+ // We're lucky to have three different MINERvA CC0pi anti-numu 2D distributions, only for LE
+ if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_antinu")) fDist = kPtPz;
else if (!name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuQE_antinu")) fDist = kQ2QEEnuQE;
else if (!name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuTrue_antinu")) fDist = kQ2QEEnuTrue;
// Define what files to use from the dist
std::string basedir = "MINERvA/CC0pi_2D_antinu/";
std::string datafile = basedir;
std::string covfile = basedir;
std::string titles = "";
std::string distdescript = "";
std::string histname = "";
std::string xbinning = basedir;
std::string ybinning = basedir;
// N.B. fScaleFactor also needs to be set dependent on the distribution. The EnuQE and EnuTrue distributions flux integrate in the Enu dimension and flux average in the Q2 dimension
switch (fDist) {
case (kPtPz):
- datafile += "cross_sections_muonpz_muonpt_lowangleqelike_minerva_2d.csv";
- covfile += "cross_sections_muonpz_muonpt_lowangleqelike_minerva_covariance.csv";
- xbinning += "cross_sections_muonpt_lowangleqelike_minerva_intmuonpz_bins_1d.csv";
- ybinning += "cross_sections_muonpz_lowangleqelike_minerva_intmuonpt_bins_1d.csv";
- titles = "MINERvA CC0#pi #bar{#nu}_{#mu} p_{t} p_{z};p_{t} (GeV);p_{z} (GeV);d^{2}#sigma/dp_{t}dp_{z} (cm^{2}/GeV^{2}/nucleon)";
distdescript = "MINERvA_CC0pi_XSec_2Dptpz_antinu sample";
+ if (IsME) {
+ datafile += "data_release_pzpt.root";
+ titles = "MINERvA CC0#pi #bar{#nu}_{#mu} p_{z} p_{t};p_{z} (GeV);p_{t} (GeV);d^{2}#sigma/dp_{z}dp_{t} (cm^{2}/GeV^{2}/nucleon)";
+ } else {
+ datafile += "cross_sections_muonpz_muonpt_lowangleqelike_minerva_2d.csv";
+ covfile += "cross_sections_muonpz_muonpt_lowangleqelike_minerva_covariance.csv";
+ xbinning += "cross_sections_muonpt_lowangleqelike_minerva_intmuonpz_bins_1d.csv";
+ ybinning += "cross_sections_muonpz_lowangleqelike_minerva_intmuonpt_bins_1d.csv";
+ titles = "MINERvA CC0#pi #bar{#nu}_{#mu} p_{t} p_{z};p_{t} (GeV);p_{z} (GeV);d^{2}#sigma/dp_{t}dp_{z} (cm^{2}/GeV^{2}/nucleon)";
+ }
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux();
break;
case (kQ2QEEnuQE):
datafile += "cross_sections_enuqe_qsq_lowangleqelike_minerva_2d.csv";
covfile += "cross_sections_enuqe_qsq_lowangleqelike_minerva_covariance.csv";
xbinning += "cross_sections_qsq_lowangleqelike_minerva_intenuqe_bins_1d.csv";
ybinning += "cross_sections_enuqe_lowangleqelike_minerva_intqsq_bins_1d.csv";
titles = "MINERvA CC0#pi #bar{#nu}_{#mu} Q^{2}_{QE} E^{#nu}_{QE};Q^{2}_{QE} (GeV);E^{#nu}_{QE} (GeV);d#sigma/dQ^{2}_{QE} (cm^{2}/GeV^{2}/nucleon)";
distdescript = "MINERvA_CC0pi_XSec_2DQ2QEEnuQE_antinu sample";
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents);
break;
case (kQ2QEEnuTrue):
datafile += "cross_sections_enu_qsq_lowangleqelike_minerva_2d.csv";
covfile += "cross_sections_enu_qsq_lowangleqelike_minerva_covariance.csv";
xbinning += "cross_sections_qsq_lowangleqelike_minerva_intenuqe_bins_1d.csv";
ybinning += "cross_sections_enuqe_lowangleqelike_minerva_intqsq_bins_1d.csv";
titles = "MINERvA CC0#pi #bar{#nu}_{#mu} Q^{2}_{QE} E^{#nu}_{True};Q^{2}_{QE} (GeV);E^{#nu}_{True} (GeV);d#sigma/dQ^{2}_{QE} (cm^{2}/GeV^{2}/nucleon)";
distdescript = "MINERvA_CC0pi_XSec_2DQ2QEEnuTrue_antinu sample";
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents);
break;
default:
NUIS_ABORT("Unknown Analysis Distribution : " << name);
}
fSettings.SetTitle( GeneralUtils::ParseToStr(titles,";")[0] );
fSettings.SetXTitle( GeneralUtils::ParseToStr(titles,";")[1] );
fSettings.SetYTitle( GeneralUtils::ParseToStr(titles,";")[2] );
fSettings.SetZTitle( GeneralUtils::ParseToStr(titles,";")[3] );
// Sample overview ---------------------------------------------------
std::string descrip = distdescript + "\n"\
- "Target: CH \n" \
- "Flux: MINERvA Low Energy RHC anti-numu \n" \
- "Signal: CC-0pi \n";
+ "Target: CH \n";
+
+ if (IsME) {
+ descrip += "Flux: MINERvA Medium Energy RHC anti-numu \n" \
+ "Signal: CC-0pi \n";
+ } else {
+ descrip += "Flux: MINERvA Low Energy RHC anti-numu \n" \
+ "Signal: CC-0pi \n";
+ }
fSettings.SetDescription(descrip);
// The input ROOT file in the fSettings
fSettings.SetDataInput(FitPar::GetDataBase() + datafile);
- fSettings.SetCovarInput(FitPar::GetDataBase() + covfile);
- // Save the binning used in the sample settings
- fSettings.SetS("xbins", FitPar::GetDataBase() + xbinning);
- fSettings.SetS("ybins", FitPar::GetDataBase() + ybinning);
-
- // Sets up the data from the data file, x binning and y binning
- SetDataFromTextFile(fSettings.GetDataInput(),
- fSettings.GetS("xbins"),
- fSettings.GetS("ybins"));
-
- // The data comes in units of 1E-41
- fDataHist->Scale(1E-41);
-
- // Setup the covariance matrix
- SetCovarFromTextFile(fSettings.GetCovarInput(), fDataHist->GetNbinsX()*fDataHist->GetNbinsY());
-
- // Set the error on the data from the covariance matrix
- StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, (TH2I*)NULL, 1.E-41, true);
-
- // In NUISANCE we assume the covar scale is 1E76 (or cross-section in 1E-38)
- // For this measurement it's actually 1E41*1E41=1E82 so need to multiply 1E82/1E76=1E6
- double ScalingFactor = 1E-3*1E-3;
- (*fFullCovar) *= ScalingFactor;
- (*covar) *= 1./ScalingFactor;
- (*fDecomp) *= 1./ScalingFactor;
+
+ if (IsME) SetDataValues(fSettings.GetDataInput(), "data_xsection_with_totErr_norm");
+ else {
+ fSettings.SetCovarInput(FitPar::GetDataBase() + covfile);
+ // Save the binning used in the sample settings
+ fSettings.SetS("xbins", FitPar::GetDataBase() + xbinning);
+ fSettings.SetS("ybins", FitPar::GetDataBase() + ybinning);
+
+ // Sets up the data from the data file, x binning and y binning
+ SetDataFromTextFile(fSettings.GetDataInput(),
+ fSettings.GetS("xbins"),
+ fSettings.GetS("ybins"));
+
+ // The data comes in units of 1E-41
+ fDataHist->Scale(1E-41);
+
+ // Setup the covariance matrix
+ SetCovarFromTextFile(fSettings.GetCovarInput(), fDataHist->GetNbinsX()*fDataHist->GetNbinsY());
+
+ // Set the error on the data from the covariance matrix
+ StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, (TH2I*)NULL, 1.E-41, true);
+
+ // In NUISANCE we assume the covar scale is 1E76 (or cross-section in 1E-38)
+ // For this measurement it's actually 1E41*1E41=1E82 so need to multiply 1E82/1E76=1E6
+ double ScalingFactor = 1E-3*1E-3;
+ (*fFullCovar) *= ScalingFactor;
+ (*covar) *= 1./ScalingFactor;
+ (*fDecomp) *= 1./ScalingFactor;
+ }
};
//********************************************************************
MINERvA_CC0pi_XSec_2D_antinu::MINERvA_CC0pi_XSec_2D_antinu(nuiskey samplekey) {
-//********************************************************************
+ //********************************************************************
+
+ IsME = false;
fSettings = LoadSampleSettings(samplekey);
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("antinumu");
// Set up the data and covariance matrix
SetupDataSettings();
FinaliseSampleSettings();
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CC0pi_XSec_2D_antinu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Checking to see if there is a Muon
if (event->NumFSParticle(-13) == 0) return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(-13)->fP;
switch (fDist) {
case (kPtPz):
{
- Double_t px = Pmu.X()/1.E3;
- Double_t py = Pmu.Y()/1.E3;
- Double_t pt = sqrt(px*px+py*py);
-
- // Don't want to assume the event generators all have neutrino coming along z
- // pz is muon momentum projected onto the neutrino direction
- Double_t pz = Pmu.Vect().Dot(Pnu.Vect()*(1.0/Pnu.Vect().Mag()))/1.E3;
-
- // Set Hist Variables
- fYVar = pz;
- fXVar = pt;
- break;
+ Double_t px = Pmu.X()/1.E3;
+ Double_t py = Pmu.Y()/1.E3;
+ Double_t pt = sqrt(px*px+py*py);
+
+ // Don't want to assume the event generators all have neutrino coming along z
+ // pz is muon momentum projected onto the neutrino direction
+ Double_t pz = Pmu.Vect().Dot(Pnu.Vect()*(1.0/Pnu.Vect().Mag()))/1.E3;
+
+ if (IsME) {
+ fYVar = pt;
+ fXVar = pz;
+ } else {
+ fYVar = pz;
+ fXVar = pt;
+ }
+ break;
}
case (kQ2QEEnuQE):
{
- double Q2qeRec = FitUtils::Q2QErec(Pmu, Pnu, 30, false);
- double EnuQErec = FitUtils::EnuQErec(Pmu, Pnu, 30, false);
- fXVar = Q2qeRec;
- fYVar = EnuQErec;
- break;
+ double Q2qeRec = FitUtils::Q2QErec(Pmu, Pnu, 30, false);
+ double EnuQErec = FitUtils::EnuQErec(Pmu, Pnu, 30, false);
+ fXVar = Q2qeRec;
+ fYVar = EnuQErec;
+ break;
}
case (kQ2QEEnuTrue):
{
- double Q2qeRec = FitUtils::Q2QErec(Pmu, Pnu, 30, false);
- double EnuTrue = Pnu.E()/1.E3;
- fXVar = Q2qeRec;
- fYVar = EnuTrue;
- break;
+ double Q2qeRec = FitUtils::Q2QErec(Pmu, Pnu, 30, false);
+ double EnuTrue = Pnu.E()/1.E3;
+ fXVar = Q2qeRec;
+ fYVar = EnuTrue;
+ break;
}
}
};
//********************************************************************
bool MINERvA_CC0pi_XSec_2D_antinu::isSignal(FitEvent *event) {
//********************************************************************
+ if (IsME) {
+ return SignalDef::isCC0pi_anti_MINERvAPTPZ_ME(event, -14, EnuMin, EnuMax);
+ }
return SignalDef::isCC0pi_anti_MINERvAPTPZ(event, -14, EnuMin, EnuMax);
};
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2D_antinu.h b/src/MINERvA/MINERvA_CC0pi_XSec_2D_antinu.h
index f884ba5..378149b 100644
--- a/src/MINERvA/MINERvA_CC0pi_XSec_2D_antinu.h
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_2D_antinu.h
@@ -1,55 +1,57 @@
// 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 .
*******************************************************************************/
#ifndef MINERVA_CC0PI_XSEC_2D_ANTINU_H_SEEN
#define MINERVA_CC0PI_XSEC_2D_ANTINU_H_SEEN
#include "Measurement2D.h"
//********************************************************************
class MINERvA_CC0pi_XSec_2D_antinu : public Measurement2D {
//********************************************************************
public:
// Constructor
MINERvA_CC0pi_XSec_2D_antinu(nuiskey samplekey);
// Destructor
virtual ~MINERvA_CC0pi_XSec_2D_antinu() {};
// Required functions
bool isSignal(FitEvent *nvect);
void FillEventVariables(FitEvent *event);
protected:
// Set up settings based on distribution
void SetupDataSettings();
// The 2D distributions we have
enum Distribution {
kPtPz,
kQ2QEEnuQE,
kQ2QEEnuTrue
};
Distribution fDist;
+
+ bool IsME;
};
#endif
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx
index 5b8d40c..428dc2c 100644
--- a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx
@@ -1,183 +1,210 @@
// 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 .
*******************************************************************************/
/*
Authors: Adrian Orea (v1 2017)
Clarence Wret (v2 2018)
*/
#include "MINERvA_CC0pi_XSec_2D_nu.h"
#include "MINERvA_SignalDef.h"
//********************************************************************
void MINERvA_CC0pi_XSec_2D_nu::SetupDataSettings() {
//********************************************************************
// Define what files to use from the dist
std::string datafile = "";
std::string corrfile = "";
std::string titles = "";
std::string distdescript = "";
std::string histname = "";
- datafile = "MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root";
- corrfile = "MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root";
- titles = "MINERvA CC0#pi #nu_{#mu} p_{t} p_{z};p_{z} (GeV);p_{t} "
- "(GeV);d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)";
- distdescript = "MINERvA_CC0pi_XSec_2Dptpz_nu sample";
- histname = "pt_pl_cross_section";
+ if (IsME) {
+ std::cerr << "Medium energy implementation does not support covariance" << std::endl;
+ std::cerr << "If you want to run anyway, comment me out" << std::endl;
+ std::cerr << "I'm in " << __FILE__ << ":" << __LINE__ << std::endl;
+ throw;
+ }
+
+ if (!IsME) {
+ datafile = "MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root";
+ corrfile = "MINERvA/CC0pi_2D/cov_ptpl_2D_qelike.root";
+ titles = "MINERvA CC0#pi #nu_{#mu} p_{t} p_{z};p_{z} (GeV);p_{t} "
+ "(GeV);d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)";
+ distdescript = "MINERvA_CC0pi_XSec_2Dptpz_nu sample";
+ histname = "pt_pl_cross_section";
+ } else {
+ datafile = "MINERvA/CC0pi_2D/minerva_cc0pi_ME_ptpz.root";
+ titles = "MINERvA CC0#pi #nu_{#mu} p_{t} p_{z};p_{z} (GeV);p_{t} "
+ "(GeV);d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)";
+ distdescript = "MINERvA_CC0pi_XSec_2Dptpz_ME_nu sample";
+ histname = "minerva_cc0pi_ME_ptpz_data";
+ }
fSettings.SetTitle(GeneralUtils::ParseToStr(titles, ";")[0]);
fSettings.SetXTitle(GeneralUtils::ParseToStr(titles, ";")[1]);
fSettings.SetYTitle(GeneralUtils::ParseToStr(titles, ";")[2]);
fSettings.SetZTitle(GeneralUtils::ParseToStr(titles, ";")[3]);
+ std::string descrip;
// Sample overview ---------------------------------------------------
- std::string descrip = distdescript + "\n"
- "Target: CH \n"
- "Flux: MINERvA Low Energy FHC numu \n"
- "Signal: CC-0pi \n";
+ if (IsME) {
+ descrip = distdescript + "\n"
+ "Target: CH \n"
+ "Flux: MINERvA Medium Energy FHC numu \n"
+ "Signal: CC-0pi \n";
+ } else {
+ descrip = distdescript + "\n"
+ "Target: CH \n"
+ "Flux: MINERvA Low Energy FHC numu \n"
+ "Signal: CC-0pi \n";
+ }
+
fSettings.SetDescription(descrip);
// The input ROOT file
fSettings.SetDataInput(FitPar::GetDataBase() + datafile);
fSettings.SetCovarInput(FitPar::GetDataBase() + corrfile);
// Set the data file
SetDataValues(fSettings.GetDataInput(), histname);
}
//********************************************************************
MINERvA_CC0pi_XSec_2D_nu::MINERvA_CC0pi_XSec_2D_nu(nuiskey samplekey) {
//********************************************************************
+ IsME = false;
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
SetupDataSettings();
FinaliseSampleSettings();
fScaleFactor =
- (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
- this->TotalIntegratedFlux();
-
- TMatrixDSym *tempmat = StatUtils::GetCovarFromRootFile(
- fSettings.GetCovarInput(), "TotalCovariance");
- fFullCovar = tempmat;
- // Decomposition is stable for entries that aren't E-xxx
- double ScalingFactor = 1E38 * 1E38;
- (*fFullCovar) *= ScalingFactor;
-
- // Just check that the data error and covariance are the same
- for (int i = 0; i < fFullCovar->GetNrows(); ++i) {
- for (int j = 0; j < fFullCovar->GetNcols(); ++j) {
- // Get the global bin
- int xbin1, ybin1, zbin1;
- fDataHist->GetBinXYZ(i, xbin1, ybin1, zbin1);
- double xlo1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1);
- double xhi1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1 + 1);
- double ylo1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1);
- double yhi1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1 + 1);
- if (xlo1 < fDataHist->GetXaxis()->GetBinLowEdge(1) ||
- ylo1 < fDataHist->GetYaxis()->GetBinLowEdge(1) ||
- xhi1 > fDataHist->GetXaxis()->GetBinLowEdge(
- fDataHist->GetXaxis()->GetNbins() + 1) ||
- yhi1 > fDataHist->GetYaxis()->GetBinLowEdge(
- fDataHist->GetYaxis()->GetNbins() + 1))
- continue;
- double data_error = fDataHist->GetBinError(xbin1, ybin1);
- double cov_error = sqrt((*fFullCovar)(i, i) / ScalingFactor);
- if (fabs(data_error - cov_error) > 1E-5) {
- std::cerr << "Error on data is different to that of covariance"
- << std::endl;
- NUIS_ERR(FTL, "Data error: " << data_error);
- NUIS_ERR(FTL, "Cov error: " << cov_error);
- NUIS_ERR(FTL, "Data/Cov: " << data_error / cov_error);
- NUIS_ERR(FTL, "Data-Cov: " << data_error - cov_error);
- NUIS_ERR(FTL, "For x: " << xlo1 << "-" << xhi1);
- NUIS_ABORT("For y: " << ylo1 << "-" << yhi1);
+ (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
+ this->TotalIntegratedFlux();
+
+ if (!IsME) {
+ TMatrixDSym *tempmat = StatUtils::GetCovarFromRootFile(
+ fSettings.GetCovarInput(), "TotalCovariance");
+ fFullCovar = tempmat;
+ // Decomposition is stable for entries that aren't E-xxx
+ double ScalingFactor = 1E38 * 1E38;
+ (*fFullCovar) *= ScalingFactor;
+
+ // Just check that the data error and covariance are the same
+ for (int i = 0; i < fFullCovar->GetNrows(); ++i) {
+ for (int j = 0; j < fFullCovar->GetNcols(); ++j) {
+ // Get the global bin
+ int xbin1, ybin1, zbin1;
+ fDataHist->GetBinXYZ(i, xbin1, ybin1, zbin1);
+ double xlo1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1);
+ double xhi1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1 + 1);
+ double ylo1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1);
+ double yhi1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1 + 1);
+ if (xlo1 < fDataHist->GetXaxis()->GetBinLowEdge(1) ||
+ ylo1 < fDataHist->GetYaxis()->GetBinLowEdge(1) ||
+ xhi1 > fDataHist->GetXaxis()->GetBinLowEdge(
+ fDataHist->GetXaxis()->GetNbins() + 1) ||
+ yhi1 > fDataHist->GetYaxis()->GetBinLowEdge(
+ fDataHist->GetYaxis()->GetNbins() + 1))
+ continue;
+ double data_error = fDataHist->GetBinError(xbin1, ybin1);
+ double cov_error = sqrt((*fFullCovar)(i, i) / ScalingFactor);
+ if (fabs(data_error - cov_error) > 1E-5) {
+ std::cerr << "Error on data is different to that of covariance"
+ << std::endl;
+ NUIS_ERR(FTL, "Data error: " << data_error);
+ NUIS_ERR(FTL, "Cov error: " << cov_error);
+ NUIS_ERR(FTL, "Data/Cov: " << data_error / cov_error);
+ NUIS_ERR(FTL, "Data-Cov: " << data_error - cov_error);
+ NUIS_ERR(FTL, "For x: " << xlo1 << "-" << xhi1);
+ NUIS_ABORT("For y: " << ylo1 << "-" << yhi1);
+ }
}
}
- }
- // Now can make the inverted covariance
- covar = StatUtils::GetInvert(fFullCovar);
- fDecomp = StatUtils::GetDecomp(fFullCovar);
+ // Now can make the inverted covariance
+ covar = StatUtils::GetInvert(fFullCovar);
+ fDecomp = StatUtils::GetDecomp(fFullCovar);
- // Now scale back
- (*fFullCovar) *= 1.0 / ScalingFactor;
- (*fDecomp) *= ScalingFactor;
- //Don't scale this back as GetLikelihood expects it to look like this
- // (*covar) *= ScalingFactor;
+ // Now scale back
+ (*fFullCovar) *= 1.0 / ScalingFactor;
+ (*fDecomp) *= ScalingFactor;
+ //Don't scale this back as GetLikelihood expects it to look like this
+ // (*covar) *= ScalingFactor;
- fMapHist = new TH2I("MINERvA_CC0pi_XSec_2D_nu_maphist", "",
- fDataHist->GetNbinsX(), 0, fDataHist->GetNbinsX(),
- fDataHist->GetNbinsY(), 0, fDataHist->GetNbinsY());
+ fMapHist = new TH2I("MINERvA_CC0pi_XSec_2D_nu_maphist", "",
+ fDataHist->GetNbinsX(), 0, fDataHist->GetNbinsX(),
+ fDataHist->GetNbinsY(), 0, fDataHist->GetNbinsY());
- int nbinsx = fDataHist->GetNbinsX();
- int nbinsy = fDataHist->GetNbinsY();
- Int_t Nbins = nbinsx * nbinsy;
+ int nbinsx = fDataHist->GetNbinsX();
+ int nbinsy = fDataHist->GetNbinsY();
+ Int_t Nbins = nbinsx * nbinsy;
- // Loop over the covariance matrix bins
- for (int i = 0; i < Nbins; ++i) {
- int xbin = (i % nbinsx) + 1;
- int ybin = (i / nbinsx) + 1;
+ // Loop over the covariance matrix bins
+ for (int i = 0; i < Nbins; ++i) {
+ int xbin = (i % nbinsx) + 1;
+ int ybin = (i / nbinsx) + 1;
- fMapHist->SetBinContent(xbin, ybin, i+1);
- }
+ fMapHist->SetBinContent(xbin, ybin, i+1);
+ }
+ }
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CC0pi_XSec_2D_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Checking to see if there is a Muon
if (event->NumFSParticle(13) == 0)
return;
// Get the muon kinematics
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
Double_t px = Pmu.X() / 1000;
Double_t py = Pmu.Y() / 1000;
Double_t pt = sqrt(px * px + py * py);
// y-axis is transverse momentum for both measurements
fYVar = pt;
// Don't want to assume the event generators all have neutrino coming along
// z pz is muon momentum projected onto the neutrino direction
Double_t pz = Pmu.Vect().Dot(Pnu.Vect() * (1.0 / Pnu.Vect().Mag())) / 1000.;
// Set Hist Variables
fXVar = pz;
};
//********************************************************************
bool MINERvA_CC0pi_XSec_2D_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCC0pi_MINERvAPTPZ(event, 14, EnuMin, EnuMax);
};
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.h b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.h
index 120fb9b..5a23030 100644
--- a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.h
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.h
@@ -1,50 +1,51 @@
// 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 .
*******************************************************************************/
#ifndef MINERVA_CC0PI_XSEC_2D_NU_H_SEEN
#define MINERVA_CC0PI_XSEC_2D_NU_H_SEEN
#include "Measurement2D.h"
//********************************************************************
class MINERvA_CC0pi_XSec_2D_nu : public Measurement2D {
//********************************************************************
public:
// Constructor
MINERvA_CC0pi_XSec_2D_nu(nuiskey samplekey);
// Destructor
virtual ~MINERvA_CC0pi_XSec_2D_nu() {};
// Required functions
bool isSignal(FitEvent *nvect);
void FillEventVariables(FitEvent *event);
protected:
// Converted covariance matrix to provide global binning method in GetLikelihood
// double GetLikelihood();
// Set up settings based on distribution
void SetupDataSettings();
+ bool IsME;
};
#endif
diff --git a/src/MINERvA/MINERvA_SignalDef.h b/src/MINERvA/MINERvA_SignalDef.h
index e0c3003..dbdc316 100644
--- a/src/MINERvA/MINERvA_SignalDef.h
+++ b/src/MINERvA/MINERvA_SignalDef.h
@@ -1,126 +1,127 @@
// 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 .
*******************************************************************************/
#ifndef MINERVA_SIGNALDEF_H_SEEN
#define MINERVA_SIGNALDEF_H_SEEN
#include "FitEvent.h"
namespace SignalDef {
// *********************************
/// MINERvA CC1pi+/- signal definition (2015 release)
/// Note: There is a 2016 release which is different to this (listed below),
/// but
/// it is CCNpi+ and has a different W cut
/// Note2: The W cut is implemented in the class implementation in MINERvA/
/// rather than here so we can draw events that don't pass the W cut (W cut is
/// 1.4 GeV)
/// Could possibly be changed for slight speed increase since less events
/// would be used
///
/// MINERvA signal is slightly different to MiniBooNE
///
/// Exactly one negative muon
/// Exactly one charged pion (both + and -); however, there is a Michel e-
/// requirement but UNCLEAR IF UNFOLDED OR NOT (so don't know if should be
/// applied)
/// Exactly 1 charged pion exits (so + and - charge), however, has Michel
/// electron requirement, so look for + only here?
/// No restriction on neutral pions or other mesons
/// MINERvA has unfolded and not unfolded muon phase space for 2015
///
/// Possible problems:
/// 1) Should there be a pi+ only cut implemented due to Michel requirement, or
/// is pi- events filled from MC?
/// 2) There is a T_pi < 350 MeV cut coming from requiring a stopping pion so
/// the
/// Michel e is seen, this is also unclear if it's unfolded so any pion is OK
///
/// Nice things:
/// Much data given: with and without muon angle cuts and with and without shape
/// only data + covariance
bool isCC1pip_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool isRestricted = false);
bool isCC1pip_MINERvA_2017(FitEvent *event, double EnuMin, double EnuMax);
// *********************************
/// MINERvA CCNpi+/- signal definition from 2016 publication
/// Different to CC1pi+/- listed above; additional has W < 1.8 GeV
///
/// Still asks for a Michel e and still unclear if this is unfolded or not
/// Says stuff like "requirement that a Michel e isolates a subsample that is
/// more nearly a pi+ prodution", yet the signal definition is both pi+ and pi-?
///
/// One negative muon
/// At least one charged pion
/// 1.5 < Enu < 10
/// No restrictions on pi0 or other mesons or baryons
///
/// Also writes number of pions (nPions) if studies on this want to be done...
bool isCCNpip_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool isRestricted = false, bool isWtrue=false);
// *********************************
/// MINERvA numubar CC 1pi- signal definition
/// Used in 2019 analysis of LE data
/// (Phys. Rev. D 100, 052008, https://doi.org/10.1103/PhysRevD.100.052008)
/// Implemented 30 April 2021 by S. Gardiner
///
/// ** Requirements on final-state particles **
/// (see Eq. (1) of the publication)
///
/// Exactly one positive muon
/// Exactly one negative pion
/// Zero additional mesons
/// Any number of nucleons
///
/// ** Kinematic limits **
/// (see Sec. V of the publication)
///
/// Muon scattering angle (theta_mu) < 25 degrees
///
/// 1.5 GeV < E_nubar < 10.0 GeV
/// (appears to be *true energy* due to unfolding)
///
/// W_exp < 1.8 GeV (experimental estimator of true hadronic invariant mass)
bool isCC1pim_MINERvA(FitEvent *event, double EnuMin, double EnuMax);
bool isCCQEnumu_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace = true);
bool isCCQEnumubar_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace = true);
bool isCCincLowRecoil_MINERvA(FitEvent *event, double EnuMin, double EnuMax);
bool isCC0pi1p_MINERvA(FitEvent *event, double enumin, double enumax);
bool isCC1pi0_MINERvA_2015(FitEvent *event, double EnuMin, double EnuMax);
bool isCC1pi0_MINERvA_2016(FitEvent *event, double EnuMin, double EnuMax);
bool isCC1pi0_MINERvA_nu(FitEvent *event, double EnuMin, double EnuMax);
bool isCC0piNp_MINERvA_STV(FitEvent *event, double EnuMin, double EnuMax);
bool isCC0pi_MINERvAPTPZ(FitEvent *event, int nuPDG, double EnuMin = 0, double EnuMax = 0);
bool isCC0pi_anti_MINERvAPTPZ(FitEvent* event, int nuPDG, double EnuMin = 0, double EnuMax = 0);
+ bool isCC0pi_anti_MINERvAPTPZ_ME(FitEvent* event, int nuPDG, double EnuMin = 0, double EnuMax = 0);
}
#endif