diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx
index 343fed5..4cffb9a 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx
@@ -1,98 +1,113 @@
#include <iomanip>
#include "T2K_SignalDef.h"
#include "T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h"
// The constructor
-T2K_CC1pip_CH_XSec_1DAdlerPhi_nu::T2K_CC1pip_CH_XSec_1DAdlerPhi_nu(nuiskey samplekey) {
+T2K_CC1pip_CH_XSec_1DAdlerPhi_nu::T2K_CC1pip_CH_XSec_1DAdlerPhi_nu(
+ nuiskey samplekey) {
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_CH_XSec_1DAdlerPhi_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2K Forward Horn Current numu \n" \
- "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n"
+ "Target: CH \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.2, costheta_pi > 0.2\n"
+ "https://arxiv.org/abs/1909.03936";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_CH_XSec_1DAdlerPhi_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#phi_{Adler} (radians)");
fSettings.SetYTitle("d#sigma/d#phi_{Adler} (cm^{2}/rad/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
- fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width");
+ fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) /
+ double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
- SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/phi_adler.rootout.root", "Phi_Adler");
- SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/phi_adler.rootout.root", "Phi_AdlerCov");
-
+ SetDataFromRootFile(GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/CH/phi_adler.rootout.root",
+ "Phi_Adler");
+ SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/CH/phi_adler.rootout.root",
+ "Phi_AdlerCov");
+
SetShapeCovar();
fDataHist->Scale(1E-38);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
-
void T2K_CC1pip_CH_XSec_1DAdlerPhi_nu::FillEventVariables(FitEvent *event) {
- if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return;
+ if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0)
+ return;
// Reconstruct the neutrino
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
- TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
+ TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
- double rEnu = FitUtils::EnuCC1piprec_T2K_eMB(Pnu, Pmu, Ppip)*1000.;
+ double rEnu = FitUtils::EnuCC1piprec_T2K_eMB(Pnu, Pmu, Ppip) * 1000.;
// Now make the reconstructed neutrino
// Has the same direction as the predicted neutrino
- TLorentzVector PnuReco(Pnu.Vect().Unit().X()*rEnu, Pnu.Vect().Unit().Y()*rEnu, Pnu.Vect().Unit().Z()*rEnu, rEnu);
+ TLorentzVector PnuReco(Pnu.Vect().Unit().X() * rEnu,
+ Pnu.Vect().Unit().Y() * rEnu,
+ Pnu.Vect().Unit().Z() * rEnu, rEnu);
// Reconstruct the initial state assuming still nucleon
- TLorentzVector Pinit(0, 0, 0, PhysConst::mass_proton*1000.);
- // Pretty much a copy of FitUtils::PhiAdler but using the reconstructed neutrino instead of true neutrino{
-
- // Get the "resonance" lorentz vector (pion proton system) reconstructed from the variables
- TLorentzVector Pres = PnuReco+Pinit-Pmu;
- // Boost the particles into the resonance rest frame so we can define the x,y,z axis
+ TLorentzVector Pinit(0, 0, 0, PhysConst::mass_proton * 1000.);
+ // Pretty much a copy of FitUtils::PhiAdler but using the reconstructed
+ // neutrino instead of true neutrino{
+
+ // Get the "resonance" lorentz vector (pion proton system) reconstructed from
+ // the variables
+ TLorentzVector Pres = PnuReco + Pinit - Pmu;
+ // Boost the particles into the resonance rest frame so we can define the
+ // x,y,z axis
PnuReco.Boost(Pres.BoostVector());
Pmu.Boost(-Pres.BoostVector());
Ppip.Boost(-Pres.BoostVector());
// The z-axis is defined as the axis of three-momentum transfer, \vec{k}
// Also unit normalise the axis
- TVector3 zAxis = (PnuReco.Vect()-Pmu.Vect())*(1.0/((PnuReco.Vect()-Pmu.Vect()).Mag()));
+ TVector3 zAxis = (PnuReco.Vect() - Pmu.Vect()) *
+ (1.0 / ((PnuReco.Vect() - Pmu.Vect()).Mag()));
- // The y-axis is then defined perpendicular to z and muon momentum in the resonance frame
+ // The y-axis is then defined perpendicular to z and muon momentum in the
+ // resonance frame
TVector3 yAxis = PnuReco.Vect().Cross(Pmu.Vect());
- yAxis *= 1.0/double(yAxis.Mag());
+ yAxis *= 1.0 / double(yAxis.Mag());
// And the x-axis is then simply perpendicular to z and x
TVector3 xAxis = yAxis.Cross(zAxis);
- xAxis *= 1.0/double(xAxis.Mag());
+ xAxis *= 1.0 / double(xAxis.Mag());
double x = Ppip.Vect().Dot(xAxis);
double y = Ppip.Vect().Dot(yAxis);
- //double z = Ppi.Vect().Dot(zAxis);
+ // double z = Ppi.Vect().Dot(zAxis);
double newphi = atan2(y, x);
// Convert negative angles to positive
- //if (newphi < 0.0) newphi += 360.0;
+ // if (newphi < 0.0) newphi += 360.0;
fXVar = newphi;
return;
};
//********************************************************************
bool T2K_CC1pip_CH_XSec_1DAdlerPhi_nu::isSignal(FitEvent *event) {
-//********************************************************************
- return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false);
+ //********************************************************************
+ return SignalDef::isCC1pip_T2K_arxiv1909_03936(
+ event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionHighEff);
}
-
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx
index d21785a..999fd7b 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx
@@ -1,86 +1,101 @@
#include <iomanip>
#include "T2K_SignalDef.h"
#include "T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h"
// The constructor
-T2K_CC1pip_CH_XSec_1DCosThAdler_nu::T2K_CC1pip_CH_XSec_1DCosThAdler_nu(nuiskey samplekey) {
+T2K_CC1pip_CH_XSec_1DCosThAdler_nu::T2K_CC1pip_CH_XSec_1DCosThAdler_nu(
+ nuiskey samplekey) {
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_CH_XSec_1DCosThAdler_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2K Forward Horn Current numu \n" \
- "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n"
+ "Target: CH \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.2, costheta_pi > 0.2\n"
+ "https://arxiv.org/abs/1909.03936";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_CH_XSec_1DCosThAdler_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("cos#theta_{Adler} (radians)");
fSettings.SetYTitle("d#sigma/dcos#theta_{Adler} (cm^{2}/1/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
- fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width");
+ fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) /
+ double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
- SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/theta_adler.rootout.root", "Theta_Adler");
- SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/theta_adler.rootout.root", "Theta_AdlerCov");
-
+ SetDataFromRootFile(GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/CH/theta_adler.rootout.root",
+ "Theta_Adler");
+ SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/CH/theta_adler.rootout.root",
+ "Theta_AdlerCov");
+
SetShapeCovar();
fDataHist->Scale(1E-38);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
-
void T2K_CC1pip_CH_XSec_1DCosThAdler_nu::FillEventVariables(FitEvent *event) {
- if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return;
+ if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0)
+ return;
- // Essentially the same as FitUtils::CosThAdler but uses the reconsturcted neutrino instead of the true neutrino
- // Reconstruct the neutrino
+ // Essentially the same as FitUtils::CosThAdler but uses the reconsturcted
+ // neutrino instead of the true neutrino Reconstruct the neutrino
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
- TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
+ TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
- double rEnu = FitUtils::EnuCC1piprec_T2K_eMB(Pnu, Pmu, Ppip)*1000.;
+ double rEnu = FitUtils::EnuCC1piprec_T2K_eMB(Pnu, Pmu, Ppip) * 1000.;
// Now make the reconstructed neutrino
// Has the same direction as the predicted neutrino
- TLorentzVector PnuReco(Pnu.Vect().Unit().X()*rEnu, Pnu.Vect().Unit().Y()*rEnu, Pnu.Vect().Unit().Z()*rEnu, rEnu);
+ TLorentzVector PnuReco(Pnu.Vect().Unit().X() * rEnu,
+ Pnu.Vect().Unit().Y() * rEnu,
+ Pnu.Vect().Unit().Z() * rEnu, rEnu);
// Reconstruct the initial state assuming still nucleon
- TLorentzVector Pinit(0, 0, 0, PhysConst::mass_proton*1000.);
- // Pretty much a copy of FitUtils::PhiAdler but using the reconstructed neutrino instead of true neutrino{
-
- // Get the "resonance" lorentz vector (pion proton system) reconstructed from the variables
- TLorentzVector Pres = PnuReco+Pinit-Pmu;
- // Boost the particles into the resonance rest frame so we can define the x,y,z axis
+ TLorentzVector Pinit(0, 0, 0, PhysConst::mass_proton * 1000.);
+ // Pretty much a copy of FitUtils::PhiAdler but using the reconstructed
+ // neutrino instead of true neutrino{
+
+ // Get the "resonance" lorentz vector (pion proton system) reconstructed from
+ // the variables
+ TLorentzVector Pres = PnuReco + Pinit - Pmu;
+ // Boost the particles into the resonance rest frame so we can define the
+ // x,y,z axis
PnuReco.Boost(Pres.BoostVector());
Pmu.Boost(-Pres.BoostVector());
Ppip.Boost(-Pres.BoostVector());
// The z-axis is defined as the axis of three-momentum transfer, \vec{k}
// Also unit normalise the axis
- TVector3 zAxis = (PnuReco.Vect()-Pmu.Vect())*(1.0/((PnuReco.Vect()-Pmu.Vect()).Mag()));
+ TVector3 zAxis = (PnuReco.Vect() - Pmu.Vect()) *
+ (1.0 / ((PnuReco.Vect() - Pmu.Vect()).Mag()));
- // Then the angle between the pion in the "resonance" rest-frame and the z-axis is the theta Adler angle
+ // Then the angle between the pion in the "resonance" rest-frame and the
+ // z-axis is the theta Adler angle
double costhAdler = cos(Ppip.Vect().Angle(zAxis));
fXVar = costhAdler;
return;
};
//********************************************************************
bool T2K_CC1pip_CH_XSec_1DCosThAdler_nu::isSignal(FitEvent *event) {
-//********************************************************************
- return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false);
+ //********************************************************************
+ return SignalDef::isCC1pip_T2K_arxiv1909_03936(
+ event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionHighEff);
}
-
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx
index d3ab741..0c16689 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx
@@ -1,118 +1,99 @@
#include <iomanip>
#include "T2K_SignalDef.h"
#include "T2K_CC1pip_CH_XSec_1DQ2_nu.h"
// The constructor
T2K_CC1pip_CH_XSec_1DQ2_nu::T2K_CC1pip_CH_XSec_1DQ2_nu(nuiskey samplekey) {
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_CH_XSec_1DQ2_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2K Forward Horn Current numu \n" \
- "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n"
+ "Target: CH \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.2, costheta_pi > 0.2\n"
+ "https://arxiv.org/abs/1909.03936";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_CH_XSec_1DQ2_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2} (GeV/c)^{2}");
fSettings.SetYTitle("d#sigma/dQ^{2} (cm^{2}/GeV^{2}/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Q2.rootout.root", "Q2");
SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Q2.rootout.root", "Q2Cov");
SetShapeCovar();
fDataHist->Scale(1E-38);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void T2K_CC1pip_CH_XSec_1DQ2_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(13) == 0 ||
event->NumFSParticle(211) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double q2 = -999;
//switch(fT2KSampleType) {
// First int refers to how we reconstruct Enu
// 0 uses true neutrino energy (not used here but common for other analyses when they unfold to true Enu from reconstructed Enu)
// 1 uses "extended MiniBooNE" method
// 2 uses "MiniBooNE reconstructed" method
// 3 uses Delta resonance mass for reconstruction
//
// The last bool refers to if pion directional information was used
//
// Use MiniBooNE reconstructed Enu; uses Michel tag so no pion direction information
//case kMB:
//q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 2, false);
//break;
// Use Extended MiniBooNE reconstructed Enu
// Needs pion information to reconstruct so bool is true (did not include Michel e tag)
//case keMB:
q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 1, true);
//break;
// Use Delta resonance reconstructed Enu
// Uses Michel electron so don't have pion directional information
//case kDelta:
//q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 3, false);
//break;
//}
fXVar = q2;
return;
};
//********************************************************************
bool T2K_CC1pip_CH_XSec_1DQ2_nu::isSignal(FitEvent *event) {
//********************************************************************
-// Warning: The CH analysis has different signal definition to the H2O analysis!
-// Often to do with the Michel tag
-
- //switch(fT2KSampleType) {
- // Using MiniBooNE formula for Enu reconstruction on the Q2 variable
- // Does have Michel e tag, set bool to true!
- //case kMB:
- //return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true);
- //break;
- // Using extended MiniBooNE formula for Enu reconstruction on the Q2 variable
- // Does not have Michel e tag because we need directional information to reconstruct Q2
- //case keMB:
- return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false);
- //break;
- // Using Delta resonance for Enu reconstruction on the Q2 variable
- // Does have Michel e tag, bool to true
- //case kDelta:
- //return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true);
- //break;
- //}
-
- // Default to return false
- //return false;
+ return SignalDef::isCC1pip_T2K_arxiv1909_03936(
+ event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionHighEff);
}
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx
index 898329c..64a5dc6 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx
@@ -1,90 +1,70 @@
#include <iomanip>
#include "T2K_SignalDef.h"
#include "T2K_CC1pip_CH_XSec_1Dppi_nu.h"
//********************************************************************
T2K_CC1pip_CH_XSec_1Dppi_nu::T2K_CC1pip_CH_XSec_1Dppi_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_CH_XSec_1Dppi_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2K Forward Horn Current numu \n" \
- "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n"
+ "Target: CH \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.2, costheta_pi > 0.2\n"
+ "https://arxiv.org/abs/1909.03936";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_CH_XSec_1Dppi_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("p_{#pi} (GeV/c)");
fSettings.SetYTitle("d#sigma/dp_{#pi} (cm^{2}/(GeV/c)/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/MomentumPion.rootout.root", "Momentum_pion");
SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/MomentumPion.rootout.root", "Momentum_pionCov");
SetShapeCovar();
fDataHist->Scale(1E-38);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void T2K_CC1pip_CH_XSec_1Dppi_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(13) == 0 ||
event->NumFSParticle(211) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double ppip = FitUtils::p(Ppip);
fXVar = ppip;
return;
};
//********************************************************************
bool T2K_CC1pip_CH_XSec_1Dppi_nu::isSignal(FitEvent *event) {
//********************************************************************
-// This distribution uses a somewhat different signal definition so might as well implement it separately here
-
- if (!SignalDef::isCC1pi(event, 14, 211, EnuMin, EnuMax)) return false;
-
- TLorentzVector Pnu = event->GetHMISParticle(14)->fP;
- TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
- TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
-
- // If this event passes the criteria on particle counting, enforce the T2K
- // ND280 phase space constraints
- // Will be different if Michel tag sample is included or not
- // Essentially, if there's a Michel tag we don't cut on the pion variables
-
- double p_mu = FitUtils::p(Pmu) * 1000;
- double cos_th_mu = cos(FitUtils::th(Pnu, Pmu));
- double cos_th_pi = cos(FitUtils::th(Pnu, Ppip));
-
- if (p_mu <= 200 || cos_th_mu <= 0.2 || cos_th_pi <= 0.2) {
- return false;
- } else {
- return true;
- }
-
- return false;
+ return SignalDef::isCC1pip_T2K_arxiv1909_03936(
+ event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionFwdHighMom);
}
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx
index 77d2048..fe6c19a 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx
@@ -1,65 +1,71 @@
#include <iomanip>
-#include "T2K_SignalDef.h"
#include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h"
+#include "T2K_SignalDef.h"
// The constructor
-T2K_CC1pip_CH_XSec_1Dthmupi_nu::T2K_CC1pip_CH_XSec_1Dthmupi_nu(nuiskey samplekey) {
+T2K_CC1pip_CH_XSec_1Dthmupi_nu::T2K_CC1pip_CH_XSec_1Dthmupi_nu(
+ nuiskey samplekey) {
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_CH_XSec_1Dthmupi_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2K Forward Horn Current numu \n" \
- "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n"
+ "Target: CH \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.2, costheta_pi > 0.2\n"
+ "https://arxiv.org/abs/1909.03936";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_CH_XSec_1Dthmupi_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#theta_{#mu,#pi} (radians)");
fSettings.SetYTitle("d#sigma/d#theta_{#mu,#pi} (cm^{2}/radians/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
- fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width");
+ fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) /
+ double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
- SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapimu.rootout.root", "Theta(pi,mu)(rads)");
- SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapimu.rootout.root", "Theta(pi,mu)(rads)Cov");
-
+ SetDataFromRootFile(GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/CH/Thetapimu.rootout.root",
+ "Theta(pi,mu)(rads)");
+ SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/CH/Thetapimu.rootout.root",
+ "Theta(pi,mu)(rads)Cov");
+
SetShapeCovar();
fDataHist->Scale(1E-38);
FinaliseMeasurement();
};
void T2K_CC1pip_CH_XSec_1Dthmupi_nu::FillEventVariables(FitEvent *event) {
- if (event->NumFSParticle(13) == 0 ||
- event->NumFSParticle(211) == 0)
+ if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0)
return;
- TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
+ TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
- TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
+ TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double thmupi = FitUtils::th(Pmu, Ppip);
fXVar = thmupi;
- //std::cout << thmupi << std::endl;
+ // std::cout << thmupi << std::endl;
return;
};
//********************************************************************
bool T2K_CC1pip_CH_XSec_1Dthmupi_nu::isSignal(FitEvent *event) {
-//********************************************************************
-// This distribution uses a somewhat different signal definition so might as well implement it separately here
- return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false);
+ //********************************************************************
+ return SignalDef::isCC1pip_T2K_arxiv1909_03936(
+ event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionHighEff);
}
-
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx
index 7473635..27b25cb 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx
@@ -1,83 +1,66 @@
#include <iomanip>
-#include "T2K_SignalDef.h"
#include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h"
+#include "T2K_SignalDef.h"
// The constructor
T2K_CC1pip_CH_XSec_1Dthpi_nu::T2K_CC1pip_CH_XSec_1Dthpi_nu(nuiskey samplekey) {
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_CH_XSec_1Dthpi_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2K Forward Horn Current numu \n" \
- "Signal: Any event with 1 muon -, 1 pion +, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n"
+ "Target: CH \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.2, costheta_pi > 0\n"
+ "https://arxiv.org/abs/1909.03936";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_CH_XSec_1Dthpi_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#theta_{#pi} (radians)");
fSettings.SetYTitle("d#sigma/d#theta_{#pi} (cm^{2}/radians/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
- fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width");
+ fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) /
+ double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
- SetDataFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapion.rootout.root", "Theta_pion");
- SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Thetapion.rootout.root", "Theta_pionCov");
-
+ SetDataFromRootFile(GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/CH/Thetapion.rootout.root",
+ "Theta_pion");
+ SetCovarFromRootFile(GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/CH/Thetapion.rootout.root",
+ "Theta_pionCov");
+
SetShapeCovar();
fDataHist->Scale(1E-38);
FinaliseMeasurement();
};
-
void T2K_CC1pip_CH_XSec_1Dthpi_nu::FillEventVariables(FitEvent *event) {
- if (event->NumFSParticle(211) == 0) return;
+ if (event->NumFSParticle(211) == 0)
+ return;
- TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
+ TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
double thpi = FitUtils::th(Pnu, Ppip);
fXVar = thpi;
};
//********************************************************************
bool T2K_CC1pip_CH_XSec_1Dthpi_nu::isSignal(FitEvent *event) {
-//********************************************************************
-// This distribution uses a somewhat different signal definition so might as well implement it separately here
-
- if (!SignalDef::isCC1pi(event, 14, 211, EnuMin, EnuMax)) return false;
-
- TLorentzVector Pnu = event->GetHMISParticle(14)->fP;
- TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
- TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
-
- // If this event passes the criteria on particle counting, enforce the T2K
- // ND280 phase space constraints
- // Will be different if Michel tag sample is included or not
- // Essentially, if there's a Michel tag we don't cut on the pion variables
-
- double p_mu = FitUtils::p(Pmu) * 1000;
- double p_pi = FitUtils::p(Ppip) * 1000;
- double cos_th_mu = cos(FitUtils::th(Pnu, Pmu));
- double cos_th_pi = cos(FitUtils::th(Pnu, Ppip));
-
- if (p_mu <= 200 || cos_th_mu <= 0.2 || cos_th_pi <= 0.0 || p_pi <= 200) {
- return false;
- } else {
- return true;
- }
-
- return false;
+ //********************************************************************
+ return SignalDef::isCC1pip_T2K_arxiv1909_03936(
+ event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionFwdHighMom);
}
-
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx
index 9b442f7..c66ddd0 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx
@@ -1,173 +1,179 @@
#include "T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h"
+#include "T2K_SignalDef.h"
//********************************************************************
-T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::T2K_CC1pip_CH_XSec_2Dpmucosmu_nu(nuiskey samplekey) {
-//********************************************************************
+T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::T2K_CC1pip_CH_XSec_2Dpmucosmu_nu(
+ nuiskey samplekey) {
+ //********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_CH_XSec_2Dpmucosmu_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2k Forward Horn Current nue + nuebar \n" \
- "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_CH_XSec_nu sample. \n"
+ "Target: CH \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, costheta_mu > 0"
+ "https://arxiv.org/abs/1909.03936";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle(" ");
- fSettings.SetYTitle("d^{2}#sigma/dp_{#mu}dcos#theta_{#mu} (cm^{2}/(GeV/c)/nucleon)");
+ fSettings.SetYTitle(
+ "d^{2}#sigma/dp_{#mu}dcos#theta_{#mu} (cm^{2}/(GeV/c)/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K_CC1pip_CH_XSec_2Dpmucosmu_nu");
- fSettings.SetDataInput( GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/PmuThetamu.root" );
- fSettings.SetCovarInput( GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/PmuThetamu.root" );
+ fSettings.SetDataInput(GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/CH/PmuThetamu.root");
+ fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/CH/PmuThetamu.root");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
- fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width");
+ fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) /
+ double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
- //SetDataValues( fSettings.GetDataInput() );
- //SetCovarMatrix( fSettings.GetCovarInput() );
+ // SetDataValues( fSettings.GetDataInput() );
+ // SetCovarMatrix( fSettings.GetCovarInput() );
SetHistograms();
- fFullCovar = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(), "Covariance_pmu_thetamu");
- covar = StatUtils::GetInvert(fFullCovar);
- fDecomp = StatUtils::GetDecomp(fFullCovar);
+ fFullCovar = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(),
+ "Covariance_pmu_thetamu");
+ covar = StatUtils::GetInvert(fFullCovar);
+ fDecomp = StatUtils::GetDecomp(fFullCovar);
SetShapeCovar();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::SetHistograms() {
TFile *data = new TFile(fSettings.GetDataInput().c_str(), "open");
- //std::string dataname = fSettings.Get
+ // std::string dataname = fSettings.Get
std::string dataname = "p_mu_theta_mu";
// Number of slices we have
const int nslices = 4;
int nbins = 0;
for (int i = 0; i < nslices; ++i) {
- TH1D* slice = (TH1D*)data->Get(Form("%s_%i", dataname.c_str(), i));
- slice = (TH1D*)slice->Clone((fName+Form("_data_slice%i", i)).c_str());
+ TH1D *slice = (TH1D *)data->Get(Form("%s_%i", dataname.c_str(), i));
+ slice = (TH1D *)slice->Clone((fName + Form("_data_slice%i", i)).c_str());
slice->Scale(1E-38);
slice->GetXaxis()->SetTitle(fSettings.GetS("xtitle").c_str());
slice->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str());
fDataHist_Slices.push_back(slice);
- fMCHist_Slices.push_back((TH1D*)slice->Clone((fName+Form("_mc_slice%i", i)).c_str()));
- SetAutoProcessTH1(fDataHist_Slices[i],kCMD_Write);
+ fMCHist_Slices.push_back(
+ (TH1D *)slice->Clone((fName + Form("_mc_slice%i", i)).c_str()));
+ SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write);
SetAutoProcessTH1(fMCHist_Slices[i]);
fMCHist_Slices[i]->Reset();
fMCHist_Slices[i]->SetLineColor(kRed);
nbins += slice->GetXaxis()->GetNbins();
}
fDataHist = new TH1D(dataname.c_str(), dataname.c_str(), nbins, 0, nbins);
fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());
int bincount = 1;
for (int i = 0; i < nslices; ++i) {
for (int j = 0; j < fDataHist_Slices[i]->GetXaxis()->GetNbins(); ++j) {
- fDataHist->SetBinContent(bincount, fDataHist_Slices[i]->GetBinContent(j+1));
- fDataHist->SetBinError(bincount, fDataHist_Slices[i]->GetBinError(j+1));
+ fDataHist->SetBinContent(bincount,
+ fDataHist_Slices[i]->GetBinContent(j + 1));
+ fDataHist->SetBinError(bincount, fDataHist_Slices[i]->GetBinError(j + 1));
TString title;
if (j == 0) {
title = "cos#theta_{#mu}=";
if (i == 0) {
title += "0-0.8, ";
} else if (i == 1) {
title += "0.8-0.85, ";
} else if (i == 2) {
title += "0.85-0.90, ";
} else if (i == 3) {
title += "0.90-1.00, ";
}
}
- title += Form("p_{#mu}=%.2f-%.2f", fDataHist_Slices[i]->GetBinLowEdge(j+1), fDataHist_Slices[i]->GetBinLowEdge(j+2));
+ title +=
+ Form("p_{#mu}=%.2f-%.2f", fDataHist_Slices[i]->GetBinLowEdge(j + 1),
+ fDataHist_Slices[i]->GetBinLowEdge(j + 2));
fDataHist->GetXaxis()->SetBinLabel(bincount, title);
bincount++;
}
}
fDataHist->GetXaxis()->SetTitle(fSettings.GetS("xtitle").c_str());
fDataHist->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str());
};
-
void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::FillEventVariables(FitEvent *event) {
- if (event->NumFSParticle(13) == 0) return;
+ if (event->NumFSParticle(13) == 0)
+ return;
- TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
- TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
+ TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
+ TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double pmu = FitUtils::p(Pmu);
double costhmu = cos(FitUtils::th(Pnu, Pmu));
fXVar = pmu;
fYVar = costhmu;
};
void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::FillHistograms() {
Measurement1D::FillHistograms();
if (Signal) {
FillMCSlice(fXVar, fYVar, Weight);
}
};
-void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::ConvertEventRates(){
+void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::ConvertEventRates() {
const int nslices = 4;
- for (int i = 0; i < nslices; i++){
+ for (int i = 0; i < nslices; i++) {
fMCHist_Slices[i]->GetSumw2();
}
// Do standard conversion.
Measurement1D::ConvertEventRates();
// First scale MC slices also by their width in Y and Z
fMCHist_Slices[0]->Scale(1.0 / 0.80);
fMCHist_Slices[1]->Scale(1.0 / 0.05);
fMCHist_Slices[2]->Scale(1.0 / 0.05);
fMCHist_Slices[3]->Scale(1.0 / 0.10);
// Now Convert into 1D list
fMCHist->Reset();
- int bincount = 1;
- for (int i = 0; i < nslices; i++){
- for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++){
- fMCHist->SetBinContent(bincount, fMCHist_Slices[i]->GetBinContent(j+1));
+ int bincount = 1;
+ for (int i = 0; i < nslices; i++) {
+ for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
+ fMCHist->SetBinContent(bincount, fMCHist_Slices[i]->GetBinContent(j + 1));
bincount++;
}
}
}
-void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::FillMCSlice(double pmu, double cosmu, double weight) {
+void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::FillMCSlice(double pmu, double cosmu,
+ double weight) {
// Hard code the bin edges in here
if (cosmu < 0.8) {
fMCHist_Slices[0]->Fill(pmu, weight);
} else if (cosmu > 0.8 && cosmu < 0.85) {
fMCHist_Slices[1]->Fill(pmu, weight);
} else if (cosmu > 0.85 && cosmu < 0.90) {
fMCHist_Slices[2]->Fill(pmu, weight);
} else if (cosmu > 0.90 && cosmu < 1.00) {
fMCHist_Slices[3]->Fill(pmu, weight);
}
};
//********************************************************************
bool T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::isSignal(FitEvent *event) {
-//********************************************************************
- if (!SignalDef::isCC1pi(event, 14, 211, EnuMin, EnuMax)) return false;
-
- TLorentzVector Pnu = event->GetHMISParticle(14)->fP;
- TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
-
- double cos_th_mu = cos(FitUtils::th(Pnu,Pmu));
- if (cos_th_mu >= 0) return true;
- return false;
+ //********************************************************************
+ return SignalDef::isCC1pip_T2K_arxiv1909_03936(event, EnuMin, EnuMax,
+ SignalDef::kMuonFwd);
};
-
diff --git a/src/T2K/T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.cxx b/src/T2K/T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.cxx
index 45c4bf1..77777d0 100644
--- a/src/T2K/T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.cxx
+++ b/src/T2K/T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.cxx
@@ -1,75 +1,77 @@
#include "T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h"
// The derived neutrino energy assuming a Delta resonance and a nucleon at rest; so only requires the outgoing muon to derive (and information on the angle between the muon and the neutrino)
// Please beware that this is NOT THE "TRUE" NEUTRINO ENERGY; IT'S A PROXY FOR THE TRUE NEUTRINO ENERGY
// Also, this is flux-integrated cross-section, not flux averaged
//********************************************************************
T2K_CC1pip_H2O_XSec_1DEnuDelta_nu::T2K_CC1pip_H2O_XSec_1DEnuDelta_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_H2O_XSec_1DEnuDelta_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2k Forward Horn Current nue + nuebar \n" \
- "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_H2O_XSec_nu sample. \n"
+ "Target: H20 \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.3, costheta_pi > 0.3\n"
+ "https://doi.org/10.1103/PhysRevD.97.012001";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_H2O_XSec_1DEnuDelta_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("E_{#nu} (GeV)");
fSettings.SetYTitle("#sigma(E_{#nu}) (cm^{2}/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetDataInput(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;EnuRec_Delta/hResultTot");
fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;EnuRec_Delta/TotalCovariance");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents);
// Plot Setup -------------------------------------------------------
SetDataFromRootFile( fSettings.GetDataInput() );
SetCovarFromRootFile( fSettings.GetCovarInput() );
ScaleCovar(1E76);
SetShapeCovar();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
// Find the muon whows kinematics we use to derive the "neutrino energy"
void T2K_CC1pip_H2O_XSec_1DEnuDelta_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Need to make sure there's a muon
if (event->NumFSParticle(13) == 0) return;
// Get the incoming neutrino
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
// Get the muon
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double Enu = FitUtils::EnuCC1piprecDelta(Pnu, Pmu);
fXVar = Enu;
return;
};
//********************************************************************
// Beware: The H2O analysis has different signal definition to the CH analysis!
bool T2K_CC1pip_H2O_XSec_1DEnuDelta_nu::isSignal(FitEvent *event) {
//********************************************************************
- return SignalDef::isCC1pip_T2K_H2O(event, EnuMin, EnuMax);
+ return SignalDef::isCC1pip_T2K_PRD97_012001(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC1pip_H2O_XSec_1DEnuMB_nu.cxx b/src/T2K/T2K_CC1pip_H2O_XSec_1DEnuMB_nu.cxx
index 1f1e902..0338dce 100644
--- a/src/T2K/T2K_CC1pip_H2O_XSec_1DEnuMB_nu.cxx
+++ b/src/T2K/T2K_CC1pip_H2O_XSec_1DEnuMB_nu.cxx
@@ -1,80 +1,82 @@
#include "T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h"
// The derived neutrino energy using the "MiniBooNE formula" (see paper)
// Essentially this is a proxy for the neutrino energy, using the outgoing pion and muon to get the reconstructed neutrino energy, assuming the struck nucleon was at rest
// Again, THIS IS NOT A "TRUE" NEUTRINO ENERGY!
//********************************************************************
T2K_CC1pip_H2O_XSec_1DEnuMB_nu::T2K_CC1pip_H2O_XSec_1DEnuMB_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_H2O_XSec_1DEnuMB_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2k Forward Horn Current nue + nuebar \n" \
- "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_H2O_XSec_nu sample. \n"
+ "Target: H20 \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.3, costheta_pi > 0.3\n"
+ "https://doi.org/10.1103/PhysRevD.97.012001";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_H2O_XSec_1DEnuMB_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("E_{#nu} (GeV)");
fSettings.SetYTitle("#sigma(E_{#nu}) (cm^{2}/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetDataInput(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;EnuRec_MB/hResultTot");
fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;EnuRec_MB/TotalCovariance");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents);
// Plot Setup -------------------------------------------------------
SetDataFromRootFile( fSettings.GetDataInput() );
SetCovarFromRootFile( fSettings.GetCovarInput() );
ScaleCovar(1E76);
SetShapeCovar();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
// Find the derived neutrino energy using the "MiniBooNE formula" (see paper)
// Essentially uses the pion and muon kinematics to derive a pseudo-neutrino energy, assuming the struck nucleon is at rest
// We also need the incoming neutrino to get the muon/neutrino and pion/neutrino angles
void T2K_CC1pip_H2O_XSec_1DEnuMB_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Need to make sure there's a muon
if (event->NumFSParticle(13) == 0) return;
// Need to make sure there's a pion
if (event->NumFSParticle(211) == 0) return;
// Get the incoming neutrino
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
// Get the muon
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
// Get the pion
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
double Enu = FitUtils::EnuCC1piprec(Pnu, Pmu, Ppip);
fXVar = Enu;
return;
};
//********************************************************************
// Beware: The H2O analysis has different signal definition to the CH analysis!
bool T2K_CC1pip_H2O_XSec_1DEnuMB_nu::isSignal(FitEvent *event) {
//********************************************************************
- return SignalDef::isCC1pip_T2K_H2O(event, EnuMin, EnuMax);
+ return SignalDef::isCC1pip_T2K_PRD97_012001(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC1pip_H2O_XSec_1Dcosmu_nu.cxx b/src/T2K/T2K_CC1pip_H2O_XSec_1Dcosmu_nu.cxx
index 8abe3e3..d8828c9 100644
--- a/src/T2K/T2K_CC1pip_H2O_XSec_1Dcosmu_nu.cxx
+++ b/src/T2K/T2K_CC1pip_H2O_XSec_1Dcosmu_nu.cxx
@@ -1,74 +1,76 @@
#include "T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h"
// The cos of the angle between the neutrino and the muon
//********************************************************************
T2K_CC1pip_H2O_XSec_1Dcosmu_nu::T2K_CC1pip_H2O_XSec_1Dcosmu_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_H2O_XSec_1Dcosmu_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2k Forward Horn Current nue + nuebar \n" \
- "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_H2O_XSec_nu sample. \n"
+ "Target: H20 \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.3, costheta_pi > 0.3\n"
+ "https://doi.org/10.1103/PhysRevD.97.012001";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_H2O_XSec_1Dcosmu_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("cos#theta_{#pi,#mu}");
fSettings.SetYTitle("d#sigma/dcos#theta_{#pi#mu} (cm^{2}/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetDataInput(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;MuCos/hResultTot");
fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;MuCos/TotalCovariance");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromRootFile( fSettings.GetDataInput() );
SetCovarFromRootFile( fSettings.GetCovarInput() );
ScaleCovar(1E76);
SetShapeCovar();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
// Find the cos theta of the angle between muon and neutrino
void T2K_CC1pip_H2O_XSec_1Dcosmu_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Need to make sure there's a muon
if (event->NumFSParticle(13) == 0) return;
// Get the incoming neutrino
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
// Get the muon
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
// Do the cos of the angle between the two
double cos_th = cos(FitUtils::th(Pnu, Pmu));
fXVar = cos_th;
return;
};
//********************************************************************
// Beware: The H2O analysis has different signal definition to the CH analysis!
bool T2K_CC1pip_H2O_XSec_1Dcosmu_nu::isSignal(FitEvent *event) {
//********************************************************************
- return SignalDef::isCC1pip_T2K_H2O(event, EnuMin, EnuMax);
+ return SignalDef::isCC1pip_T2K_PRD97_012001(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.cxx b/src/T2K/T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.cxx
index 1a50ac3..1760f7a 100644
--- a/src/T2K/T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.cxx
+++ b/src/T2K/T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.cxx
@@ -1,72 +1,74 @@
#include "T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h"
//********************************************************************
T2K_CC1pip_H2O_XSec_1Dcosmupi_nu::T2K_CC1pip_H2O_XSec_1Dcosmupi_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_H2O_XSec_1Dcosmupi_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2k Forward Horn Current nue + nuebar \n" \
- "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_H2O_XSec_nu sample. \n"
+ "Target: H20 \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.3, costheta_pi > 0.3\n"
+ "https://doi.org/10.1103/PhysRevD.97.012001";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_H2O_XSec_1Dcosmupi_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("cos#theta_{#pi,#mu}");
fSettings.SetYTitle("d#sigma/dcos#theta_{#pi#mu} (cm^{2}/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetDataInput(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;MuPiCos/hResultTot");
fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;MuPiCos/TotalCovariance");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width")*1E-38/double(fNEvents)/TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromRootFile( fSettings.GetDataInput() );
SetCovarFromRootFile( fSettings.GetCovarInput() );
ScaleCovar(1E76);
SetShapeCovar();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
// Find the cos theta of the angle between muon and pion
void T2K_CC1pip_H2O_XSec_1Dcosmupi_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Need to make sure there's a muon
if (event->NumFSParticle(13) == 0) return;
// Need to make sure there's a pion
if (event->NumFSParticle(211) == 0) return;
// Get the muon
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
// Get the pion
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
double cos_th = cos(FitUtils::th(Pmu, Ppip));
fXVar = cos_th;
return;
};
//********************************************************************
// Beware: The H2O analysis has different signal definition to the CH analysis!
bool T2K_CC1pip_H2O_XSec_1Dcosmupi_nu::isSignal(FitEvent *event) {
//********************************************************************
- return SignalDef::isCC1pip_T2K_H2O(event, EnuMin, EnuMax);
+ return SignalDef::isCC1pip_T2K_PRD97_012001(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC1pip_H2O_XSec_1Dcospi_nu.cxx b/src/T2K/T2K_CC1pip_H2O_XSec_1Dcospi_nu.cxx
index fd022ff..de63d80 100644
--- a/src/T2K/T2K_CC1pip_H2O_XSec_1Dcospi_nu.cxx
+++ b/src/T2K/T2K_CC1pip_H2O_XSec_1Dcospi_nu.cxx
@@ -1,69 +1,71 @@
#include "T2K_CC1pip_H2O_XSec_1Dcospi_nu.h"
//********************************************************************
T2K_CC1pip_H2O_XSec_1Dcospi_nu::T2K_CC1pip_H2O_XSec_1Dcospi_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_H2O_XSec_1Dcospi_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2k Forward Horn Current nue + nuebar \n" \
- "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_H2O_XSec_nu sample. \n"
+ "Target: H20 \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.3, costheta_pi > 0.3\n"
+ "https://doi.org/10.1103/PhysRevD.97.012001";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_H2O_XSec_1Dcospi_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("cos#theta_{#pi}");
fSettings.SetYTitle("d#sigma/dcos#theta_{#pi} (cm^{2}/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetDataInput(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;PosPionCos/hResultTot");
fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;PosPionCos/TotalCovariance");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromRootFile( fSettings.GetDataInput() );
SetCovarFromRootFile( fSettings.GetCovarInput() );
ScaleCovar(1E76);
SetShapeCovar();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
// Find the cos theta of the angle between pion and neutrino
void T2K_CC1pip_H2O_XSec_1Dcospi_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Need to make sure there's a pion
if (event->NumFSParticle(211) == 0) return;
// Get the incoming neutrino
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
// Get the pion
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
double cos_th = cos(FitUtils::th(Pnu, Ppip));
fXVar = cos_th;
return;
};
//********************************************************************
// Beware: The H2O analysis has different signal definition to the CH analysis!
bool T2K_CC1pip_H2O_XSec_1Dcospi_nu::isSignal(FitEvent *event) {
//********************************************************************
- return SignalDef::isCC1pip_T2K_H2O(event, EnuMin, EnuMax);
+ return SignalDef::isCC1pip_T2K_PRD97_012001(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC1pip_H2O_XSec_1Dpmu_nu.cxx b/src/T2K/T2K_CC1pip_H2O_XSec_1Dpmu_nu.cxx
index 9d6d983..48ff5ee 100644
--- a/src/T2K/T2K_CC1pip_H2O_XSec_1Dpmu_nu.cxx
+++ b/src/T2K/T2K_CC1pip_H2O_XSec_1Dpmu_nu.cxx
@@ -1,71 +1,73 @@
#include "T2K_CC1pip_H2O_XSec_1Dpmu_nu.h"
// The muon momentum
//********************************************************************
T2K_CC1pip_H2O_XSec_1Dpmu_nu::T2K_CC1pip_H2O_XSec_1Dpmu_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_H2O_XSec_1Dpmu_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2k Forward Horn Current nue + nuebar \n" \
- "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_H2O_XSec_nu sample. \n"
+ "Target: H20 \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.3, costheta_pi > 0.3\n"
+ "https://doi.org/10.1103/PhysRevD.97.012001";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_H2O_XSec_1Dpmu_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("E_{#nu} (GeV)");
fSettings.SetYTitle("#sigma(E_{#nu}) (cm^{2}/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetDataInput(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;MuMom/hResultTot");
fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;MuMom/TotalCovariance");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromRootFile( fSettings.GetDataInput() );
SetCovarFromRootFile( fSettings.GetCovarInput() );
ScaleCovar(1E76);
SetShapeCovar();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
// Find the momentum of the muon
void T2K_CC1pip_H2O_XSec_1Dpmu_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
// Need to make sure there's a muon
if (event->NumFSParticle(13) == 0) return;
// Get the muon
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double p_mu = FitUtils::p(Pmu);
fXVar = p_mu;
return;
};
//********************************************************************
// Beware: The H2O analysis has different signal definition to the CH analysis!
bool T2K_CC1pip_H2O_XSec_1Dpmu_nu::isSignal(FitEvent *event) {
//********************************************************************
- return SignalDef::isCC1pip_T2K_H2O(event, EnuMin, EnuMax);
+ return SignalDef::isCC1pip_T2K_PRD97_012001(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx b/src/T2K/T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx
index 7b33785..246c0ca 100644
--- a/src/T2K/T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx
+++ b/src/T2K/T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx
@@ -1,70 +1,78 @@
#include "T2K_CC1pip_H2O_XSec_1Dppi_nu.h"
// The momentum of the (positive) pion
//********************************************************************
T2K_CC1pip_H2O_XSec_1Dppi_nu::T2K_CC1pip_H2O_XSec_1Dppi_nu(nuiskey samplekey) {
-//********************************************************************
+ //********************************************************************
// Sample overview ---------------------------------------------------
- std::string descrip = "T2K_CC1pip_H2O_XSec_1Dppi_nu sample. \n" \
- "Target: CH \n" \
- "Flux: T2k Forward Horn Current nue + nuebar \n" \
- "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
+ std::string descrip = "T2K_CC1pip_H2O_XSec_nu sample. \n"
+ "Target: H20 \n"
+ "Flux: T2K FHC numu \n"
+ "Signal: CC1pi+, p_mu > 200 MeV, p_pi > 200 MeV\n"
+ ", costheta_mu > 0.3, costheta_pi > 0.3\n"
+ "https://doi.org/10.1103/PhysRevD.97.012001";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("T2K_CC1pip_H2O_XSec_1Dppi_nu");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("p_{#pi^{+}} (GeV/c)");
fSettings.SetYTitle("d#sigma/dp_{#pi^{+}} (cm^{2}/(GeV/c)/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
- fSettings.SetDataInput(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;PosPionMom/hResultTot");
- fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/H2O/nd280data-numu-cc1pi-xs-on-h2o-2015.root;PosPionMom/TotalCovariance");
+ fSettings.SetDataInput(
+ GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/H2O/"
+ "nd280data-numu-cc1pi-xs-on-h2o-2015.root;PosPionMom/hResultTot");
+ fSettings.SetCovarInput(
+ GeneralUtils::GetTopLevelDir() +
+ "/data/T2K/CC1pip/H2O/"
+ "nd280data-numu-cc1pi-xs-on-h2o-2015.root;PosPionMom/TotalCovariance");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
- fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width");
+ fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) /
+ double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
- SetDataFromRootFile( fSettings.GetDataInput() );
- SetCovarFromRootFile( fSettings.GetCovarInput() );
+ SetDataFromRootFile(fSettings.GetDataInput());
+ SetCovarFromRootFile(fSettings.GetCovarInput());
ScaleCovar(1E76);
SetShapeCovar();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
-
};
-
//********************************************************************
// Find the momentum of the (positively charged) pion
void T2K_CC1pip_H2O_XSec_1Dppi_nu::FillEventVariables(FitEvent *event) {
-//********************************************************************
+ //********************************************************************
// Need to make sure there's a muon
- if (event->NumFSParticle(211) == 0) return;
+ if (event->NumFSParticle(211) == 0)
+ return;
// Get the pion
- TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
+ TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
double p_pi = FitUtils::p(Ppip);
fXVar = p_pi;
return;
};
//********************************************************************
// Beware: The H2O analysis has different signal definition to the CH analysis!
bool T2K_CC1pip_H2O_XSec_1Dppi_nu::isSignal(FitEvent *event) {
-//********************************************************************
- return SignalDef::isCC1pip_T2K_H2O(event, EnuMin, EnuMax);
+ //********************************************************************
+ return SignalDef::isCC1pip_T2K_PRD97_012001(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx b/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx
index 5624011..06fbefb 100644
--- a/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx
+++ b/src/T2K/T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx
@@ -1,256 +1,257 @@
#include "T2K_CCinc_XSec_2DPcos_nu_nonuniform.h"
#include "T2K_SignalDef.h"
// ***********************************
// Implemented by Alfonso Garcia, Barcelona (now NIKHEF)
// Clarence Wret, Rochester
// (Alfonso was the T2K analyser)
// ***********************************
//********************************************************************
T2K_CCinc_XSec_2DPcos_nu_nonuniform::T2K_CCinc_XSec_2DPcos_nu_nonuniform(
nuiskey samplekey) {
//********************************************************************
fAllowedTypes += "/GENIE/NEUT";
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CCinc_XSec_2DPcos_nu_nonuniform sample. \n"
"Target: CH \n"
"Flux: T2K FHC numu \n"
- "Signal: CC-inclusive \n";
+ "Signal: CC-inclusive \n"
+ "https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.012004";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("");
// fSettings.SetXTitle("p_{#mu} (GeV)");
// fSettings.SetYTitle("cos#theta_{#mu}");
fSettings.SetYTitle("#frac{d^{2}#sigma}{dp_{#mu}dcos#theta_{#mu}} "
"[#frac{cm^{2}}{nucleon/GeV/c}]");
fSettings.SetEnuRange(0.0, 30.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K CC-inclusive p_{#mu} cos#theta_{#mu}");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * 1E-38 / fNEvents /
TotalIntegratedFlux();
// Default to using the NEUT unfolded data
UnfoldWithGENIE = false;
// Check option
if (fSettings.Found("type", "GENIE"))
UnfoldWithGENIE = true;
// Tell user what's happening
if (UnfoldWithGENIE) {
NUIS_LOG(SAM, fName << " is using GENIE unfolded data. Want NEUT? Specify "
"type=\"NEUT\" in your config file");
} else {
NUIS_LOG(SAM, fName << " is using NEUT unfolded data. Want GENIE? Specify "
"type=\"GENIE\" in your config file");
}
// Setup Histograms
SetHistograms();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
// Signal is simply a CC inclusive without any angular/momentum cuts
bool T2K_CCinc_XSec_2DPcos_nu_nonuniform::isSignal(FitEvent *event) {
return SignalDef::isCCINC(event, 14, EnuMin, EnuMax);
};
void T2K_CCinc_XSec_2DPcos_nu_nonuniform::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double pmu = Pmu.Vect().Mag() / 1000.;
double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect()));
fXVar = pmu;
fYVar = CosThetaMu;
};
// Fill up the MCSlice
void T2K_CCinc_XSec_2DPcos_nu_nonuniform::FillHistograms() {
if (Signal)
FillMCSlice(fXVar, fYVar, Weight);
}
// Modification is needed after the full reconfigure to move bins around
void T2K_CCinc_XSec_2DPcos_nu_nonuniform::ConvertEventRates() {
// Do standard conversion.
Measurement1D::ConvertEventRates();
// First scale MC slices also by their width in Y
fMCHist_Slices[0]->Scale(1.0 / 0.25);
fMCHist_Slices[1]->Scale(1.0 / 0.50);
fMCHist_Slices[2]->Scale(1.0 / 0.20);
fMCHist_Slices[3]->Scale(1.0 / 0.15);
fMCHist_Slices[4]->Scale(1.0 / 0.11);
fMCHist_Slices[5]->Scale(1.0 / 0.09);
fMCHist_Slices[6]->Scale(1.0 / 0.07);
fMCHist_Slices[7]->Scale(1.0 / 0.05);
fMCHist_Slices[8]->Scale(1.0 / 0.04);
fMCHist_Slices[9]->Scale(1.0 / 0.025);
fMCHist_Slices[10]->Scale(1.0 / 0.015);
// Now Convert into 1D list
fMCHist->Reset();
int bincount = 0;
for (int i = 0; i < nSlices; i++) {
for (int j = 0; j < fMCHist_Slices[i]->GetNbinsX(); j++) {
fMCHist->SetBinContent(bincount + 1,
fMCHist_Slices[i]->GetBinContent(j + 1));
fMCHist->SetBinError(bincount + 1, fMCHist_Slices[i]->GetBinError(j + 1));
bincount++;
}
}
};
void T2K_CCinc_XSec_2DPcos_nu_nonuniform::FillMCSlice(double x, double y,
double w) {
if (y >= -1.0 && y < -0.25)
fMCHist_Slices[0]->Fill(x, w);
else if (y >= -0.25 && y < 0.25)
fMCHist_Slices[1]->Fill(x, w);
else if (y >= 0.25 && y < 0.45)
fMCHist_Slices[2]->Fill(x, w);
else if (y >= 0.45 && y < 0.6)
fMCHist_Slices[3]->Fill(x, w);
else if (y >= 0.6 && y < 0.71)
fMCHist_Slices[4]->Fill(x, w);
else if (y >= 0.71 && y < 0.80)
fMCHist_Slices[5]->Fill(x, w);
else if (y >= 0.80 && y < 0.87)
fMCHist_Slices[6]->Fill(x, w);
else if (y >= 0.87 && y < 0.92)
fMCHist_Slices[7]->Fill(x, w);
else if (y >= 0.92 && y < 0.96)
fMCHist_Slices[8]->Fill(x, w);
else if (y >= 0.96 && y <= 0.985)
fMCHist_Slices[9]->Fill(x, w);
else if (y >= 0.985 && y <= 1.0)
fMCHist_Slices[10]->Fill(x, w);
};
void T2K_CCinc_XSec_2DPcos_nu_nonuniform::SetHistograms() {
// Read in 1D Data Histograms
TFile *fInputFile =
new TFile((FitPar::GetDataBase() +
"T2K/CCinc/nd280data-numu-cc-inc-xs-on-c-2018/histograms.root")
.c_str(),
"OPEN");
// Number of theta slices in the release
nSlices = 11;
// Data release includes unfolding with NEUT or GENIE as prior
// Choose whichever the user specifies
std::string basename;
if (UnfoldWithGENIE)
basename = "hist_xsec_data_prior_neut_cthbin";
else
basename = "hist_xsec_data_prior_neut_cthbin";
// Read in 2D Data Slices and Make MC Slices
// Count the number of bins we have in total so we can match covariance matrix
int bincount = 0;
for (int i = 0; i < nSlices; i++) {
// Get Data Histogram
// fDataHist_Slices.push_back((TH1D*)fInputFile->Get(Form("dataslice_%i",i))->Clone());
fDataHist_Slices.push_back(
(TH1D *)fInputFile->Get(Form("%s%i", basename.c_str(), i))
->Clone(
Form("T2K_CCinc_XSec_2DPcos_nu_nonuniform_slice%i_data", i)));
fDataHist_Slices[i]->SetDirectory(0);
fDataHist_Slices[i]->Scale(1E-39);
fDataHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str());
// Count up the bins
bincount += fDataHist_Slices.back()->GetNbinsX();
// Make MC Clones
fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone(
Form("T2K_CCinc_XSec_2DPcos_nu_nonuniform_Slice%i_MC", i)));
fMCHist_Slices[i]->Reset();
fMCHist_Slices[i]->SetDirectory(0);
fMCHist_Slices[i]->SetLineColor(kRed);
fMCHist_Slices[i]->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str());
SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write);
SetAutoProcessTH1(fMCHist_Slices[i]);
}
fDataHist =
new TH1D((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str(), bincount, 0, bincount);
fDataHist->SetDirectory(0);
int counter = 0;
for (int i = 0; i < nSlices; ++i) {
// Set a nice title
std::string costitle = fDataHist_Slices[i]->GetTitle();
costitle = costitle.substr(costitle.find("-> ") + 3, costitle.size());
std::string found = costitle.substr(0, costitle.find(" < "));
std::string comp = costitle.substr(costitle.find(found) + found.size() + 3,
costitle.size());
comp = comp.substr(comp.find(" < ") + 3, comp.size());
costitle = "cos#theta_{#mu}=" + found + "-" + comp;
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
fDataHist->SetBinContent(counter + 1,
fDataHist_Slices[i]->GetBinContent(j + 1));
fDataHist->SetBinError(counter + 1,
fDataHist_Slices[i]->GetBinError(j + 1));
// Set a nice axis
if (j == 0)
fDataHist->GetXaxis()->SetBinLabel(
counter + 1, Form("%s, p_{#mu}=%.1f-%.1f", costitle.c_str(),
fDataHist_Slices[i]->GetBinLowEdge(j + 1),
fDataHist_Slices[i]->GetBinLowEdge(j + 2)));
else
fDataHist->GetXaxis()->SetBinLabel(
counter + 1,
Form("p_{#mu}=%.1f-%.1f", fDataHist_Slices[i]->GetBinLowEdge(j + 1),
fDataHist_Slices[i]->GetBinLowEdge(j + 2)));
counter++;
}
}
// The correlation matrix
// Major in angular bins, minor in momentum bins: runs theta1, pmu1, pmu2,
// theta2, pmu1, pmu2, theta3, pmu1, pmu2 etc The correlation matrix
TH2D *tempcov = NULL;
if (UnfoldWithGENIE)
tempcov = (TH2D *)fInputFile->Get("covariance_matrix_genie");
else
tempcov = (TH2D *)fInputFile->Get("covariance_matrix_neut");
fFullCovar = new TMatrixDSym(bincount);
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
(*fFullCovar) *= 1E38 * 1E38;
(*covar) *= 1E-38 * 1E-38;
fInputFile->Close();
};
diff --git a/src/T2K/T2K_SignalDef.cxx b/src/T2K/T2K_SignalDef.cxx
index f6b4e7d..b407cb6 100644
--- a/src/T2K/T2K_SignalDef.cxx
+++ b/src/T2K/T2K_SignalDef.cxx
@@ -1,310 +1,343 @@
// 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/>.
-*******************************************************************************/
+ * 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 "FitUtils.h"
#include "T2K_SignalDef.h"
+#include "FitUtils.h"
namespace SignalDef {
// T2K H2O signal definition
-bool isCC1pip_T2K_H2O(FitEvent *event, double EnuMin,
- double EnuMax) {
+// https://doi.org/10.1103/PhysRevD.97.012001
+bool isCC1pip_T2K_PRD97_012001(FitEvent *event, double EnuMin, double EnuMax) {
- if (!isCC1pi(event, 14, 211, EnuMin, EnuMax)) return false;
+ if (!isCC1pi(event, 14, 211, EnuMin, EnuMax))
+ return false;
TLorentzVector Pnu = event->GetHMISParticle(14)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
double p_mu = FitUtils::p(Pmu) * 1000;
double p_pi = FitUtils::p(Ppip) * 1000;
double cos_th_mu = cos(FitUtils::th(Pnu, Pmu));
double cos_th_pi = cos(FitUtils::th(Pnu, Ppip));
if (p_mu <= 200 || p_pi <= 200 || cos_th_mu <= 0.3 || cos_th_pi <= 0.3) {
return false;
}
return true;
};
// ******************************************************
// T2K CC1pi+ CH analysis (Raquel's thesis)
// Has different phase space cuts depending on if using Michel tag or not
//
// Essentially consists of two samples: one sample which has Michel e (which we
// can't get pion direction from); this covers backwards angles quite well.
// Measurements including this sample does not have include pion kinematics cuts
// one sample which has PID in FGD and TPC
// and no Michel e. These are mostly
// forward-going so require a pion
// kinematics cut
//
// Essentially, cuts are:
// 1 negative muon
// 1 positive pion
// Any number of nucleons
// No other particles in the final state
//
-// MOVE T2K
-bool isCC1pip_T2K_CH(FitEvent *event, double EnuMin, double EnuMax,
- bool MichelElectron) {
+// https://arxiv.org/abs/1909.03936
+bool isCC1pip_T2K_arxiv1909_03936(FitEvent *event, double EnuMin,
+ double EnuMax,
+ int pscuts) {
// ******************************************************
- if (!isCC1pi(event, 14, 211, EnuMin, EnuMax)) return false;
+ if (!isCC1pi(event, 14, 211, EnuMin, EnuMax)) {
+ return false;
+ }
TLorentzVector Pnu = event->GetHMISParticle(14)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
- TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
- // If this event passes the criteria on particle counting, enforce the T2K
- // ND280 phase space constraints
- // Will be different if Michel tag sample is included or not
- // Essentially, if there's a Michel tag we don't cut on the pion variables
+ double cos_th_mu = cos(FitUtils::th(Pnu, Pmu));
+
+ if (pscuts == kMuonFwd) {
+ return (cos_th_mu > 0);
+ }
double p_mu = FitUtils::p(Pmu) * 1000;
- double p_pi = FitUtils::p(Ppip) * 1000;
- double cos_th_mu = cos(FitUtils::th(Pnu, Pmu));
- double cos_th_pi = cos(FitUtils::th(Pnu, Ppip));
- // If we're using Michel e- requirement we only care about the muon restricted
- // phase space and use full pion phase space
- if (MichelElectron) {
- // Make the cuts on the muon variables
- if (p_mu <= 200 || cos_th_mu <= 0.2) {
+ if (pscuts & kMuonHighEff) {
+ if ((cos_th_mu <= 0.2) || (p_mu <= 200)) {
return false;
- } else {
- return true;
}
+ }
- // If we aren't using Michel e- (i.e. we use directional information from
- // pion) we need to impose the phase space cuts on the muon AND the pion)
- } else {
- // Make the cuts on muon and pion variables
- if (p_mu <= 200 || p_pi <= 200 || cos_th_mu <= 0.2 || cos_th_pi <= 0.2) {
- return false;
- } else {
- return true;
- }
+ int npicuts = 0;
+ npicuts += bool(pscuts & kPionFwdHighMom);
+ npicuts += bool(pscuts & kPionHighMom);
+ npicuts += bool(pscuts & kPionHighEff);
+
+ if (npicuts != 1) {
+ NUIS_ABORT(
+ "isCC1pip_T2K_arxiv1909_03936 signal definition passed incompatible "
+ "pion phase space cuts. Should be either kMuonHighEff, or one of "
+ "kPionFwdHighMom,kPionHighMom, or kPionHighEff");
+ }
+
+ TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
+ double cos_th_pi = cos(FitUtils::th(Pnu, Ppip));
+ double p_pi = FitUtils::p(Ppip) * 1000;
+
+ if (pscuts & kPionFwdHighMom) {
+ return ((cos_th_pi > 0) && (p_pi > 200));
+ }
+
+ if (pscuts & kPionHighMom) {
+ return (p_pi > 200);
+ }
+
+ if (pscuts & kMuonHighEff) {
+ return ((cos_th_pi > 0.0) && (p_pi > 200));
}
- // Default to false; should never fire
return false;
};
bool isT2K_CC0pi(FitEvent *event, double EnuMin, double EnuMax,
- bool forwardgoing) {
+ bool forwardgoing) {
// Require a numu CC0pi event
- if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false;
+ if (!isCC0pi(event, 14, EnuMin, EnuMax))
+ return false;
TLorentzVector Pnu = event->GetHMISParticle(14)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect()));
double p_mu = Pmu.Vect().Mag();
// If we're doing a restricted phase space, Analysis II asks for:
// Cos(theta_mu) > 0.0 and p_mu > 200 MeV
if (forwardgoing) {
- if (CosThetaMu < 0.0 || p_mu < 200) return false;
+ if (CosThetaMu < 0.0 || p_mu < 200)
+ return false;
}
return true;
}
bool isT2K_CC0pi1p(FitEvent *event, double EnuMin, double EnuMax) {
// Require a numu CC0pi event
- if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false;
+ if (!isCC0pi(event, 14, EnuMin, EnuMax))
+ return false;
// Require at least one FS proton
- if (event->NumFSParticle(2212) == 0) return false;
+ if (event->NumFSParticle(2212) == 0)
+ return false;
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
- TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
-
+ TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
// Proton phase space
if (pp.Vect().Mag() < 500) {
return false;
}
- //Need exactly one proton with 500 MeV or more momentum
- std::vector<FitParticle*> protons = event->GetAllFSProton();
- int nProtonsAboveThresh=0;
- for(size_t i=0; i<protons.size(); i++){
- if(protons[i]->p()>500) nProtonsAboveThresh++;
+ // Need exactly one proton with 500 MeV or more momentum
+ std::vector<FitParticle *> protons = event->GetAllFSProton();
+ int nProtonsAboveThresh = 0;
+ for (size_t i = 0; i < protons.size(); i++) {
+ if (protons[i]->p() > 500)
+ nProtonsAboveThresh++;
}
- if(nProtonsAboveThresh!=1) return false;
+ if (nProtonsAboveThresh != 1)
+ return false;
return true;
}
-
-//CC0pi antinu in the P0D - TN328
+// CC0pi antinu in the P0D - TN328
bool isT2K_CC0piAnuP0D(FitEvent *event, double EnuMin, double EnuMax) {
// Require a anumu CC0pi event
- if (!isCC0pi(event, -14, EnuMin, EnuMax)) return false;
-
+ if (!isCC0pi(event, -14, EnuMin, EnuMax))
+ return false;
TLorentzVector pnu = event->GetHMISParticle(-14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(-13)->fP;
double Pmu = pmu.Vect().Mag();
double CosThetaMu = cos(pnu.Vect().Angle(pmu.Vect()));
// Muon phase space
- if (Pmu < 400 || Pmu > 3410) return false;
- if (Pmu < 530 && CosThetaMu<0.85) return false;
- if (Pmu < 670 && CosThetaMu<0.88) return false;
- if (Pmu < 800 && CosThetaMu<0.9) return false;
- if (Pmu < 1000 && CosThetaMu<0.91) return false;
- if (Pmu < 1380 && CosThetaMu<0.92) return false;
- if (Pmu < 2010 && CosThetaMu<0.95) return false;
-
-
+ if (Pmu < 400 || Pmu > 3410)
+ return false;
+ if (Pmu < 530 && CosThetaMu < 0.85)
+ return false;
+ if (Pmu < 670 && CosThetaMu < 0.88)
+ return false;
+ if (Pmu < 800 && CosThetaMu < 0.9)
+ return false;
+ if (Pmu < 1000 && CosThetaMu < 0.91)
+ return false;
+ if (Pmu < 1380 && CosThetaMu < 0.92)
+ return false;
+ if (Pmu < 2010 && CosThetaMu < 0.95)
+ return false;
+
return true;
}
bool isT2K_CC0piNp(FitEvent *event, double EnuMin, double EnuMax) {
// Require a numu CC0pi event
- if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false;
+ if (!isCC0pi(event, 14, EnuMin, EnuMax))
+ return false;
// Require at least one FS proton
- if (event->NumFSParticle(2212) == 0) return false;
+ if (event->NumFSParticle(2212) == 0)
+ return false;
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
- TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
+ TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
// Proton phase space
if (pp.Vect().Mag() < 500) {
return false;
}
- //Need exactly one proton with 500 MeV or more momentum
- std::vector<FitParticle*> protons = event->GetAllFSProton();
- int nProtonsAboveThresh=0;
- for(size_t i=0; i<protons.size(); i++){
- if(protons[i]->p()>500) nProtonsAboveThresh++;
+ // Need exactly one proton with 500 MeV or more momentum
+ std::vector<FitParticle *> protons = event->GetAllFSProton();
+ int nProtonsAboveThresh = 0;
+ for (size_t i = 0; i < protons.size(); i++) {
+ if (protons[i]->p() > 500)
+ nProtonsAboveThresh++;
}
- if(nProtonsAboveThresh<2) return false;
+ if (nProtonsAboveThresh < 2)
+ return false;
return true;
}
bool isT2K_CC0pi0p(FitEvent *event, double EnuMin, double EnuMax) {
// Require a numu CC0pi event
- if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false;
+ if (!isCC0pi(event, 14, EnuMin, EnuMax))
+ return false;
// Require at least one FS proton
- if (event->NumFSParticle(2212) == 0) return false;
+ if (event->NumFSParticle(2212) == 0)
+ return false;
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
- TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
-
+ TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
// Proton phase space
if (pp.Vect().Mag() > 500) {
return false;
}
return true;
}
-
bool isT2K_CC0pi_STV(FitEvent *event, double EnuMin, double EnuMax) {
// Require a numu CC0pi event
- if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false;
+ if (!isCC0pi(event, 14, EnuMin, EnuMax))
+ return false;
// Require at least one FS proton
- if (event->NumFSParticle(2212) == 0) return false;
+ if (event->NumFSParticle(2212) == 0)
+ return false;
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
- TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
+ TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
// Muon phase space
// Pmu > 250 MeV, cos(theta_mu) > -0.6 (Sweet phase space!)
if ((pmu.Vect().Mag() < 250) || cos(pnu.Vect().Angle(pmu.Vect())) < -0.6) {
return false;
}
// Proton phase space
// Pprot > 450 MeV, cos(theta_proton) > 0.4
if ((pp.Vect().Mag() < 450) || (pp.Vect().Mag() > 1E3) ||
(cos(pnu.Vect().Angle(pp.Vect())) < 0.4)) {
return false;
}
return true;
}
bool isT2K_CC0pi_ifk(FitEvent *event, double EnuMin, double EnuMax) {
// Require a numu CC0pi event
- if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false;
+ if (!isCC0pi(event, 14, EnuMin, EnuMax))
+ return false;
// Require at least one FS proton
- if (event->NumFSParticle(2212) == 0) return false;
+ if (event->NumFSParticle(2212) == 0)
+ return false;
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
- TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
+ TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
// Proton phase space
// Pprot > 450 MeV, cos(theta_proton) > 0.4
if ((pp.Vect().Mag() < 450) || (cos(pnu.Vect().Angle(pp.Vect())) < 0.4)) {
return false;
}
return true;
}
bool isT2K_CC0pi_1bin(FitEvent *event, double EnuMin, double EnuMax) {
// Require a numu CC0pi event
- if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false;
+ if (!isCC0pi(event, 14, EnuMin, EnuMax))
+ return false;
// Require at least one FS proton
- if (event->NumFSParticle(2212) == 0) return false;
+ if (event->NumFSParticle(2212) == 0)
+ return false;
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
- TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
+ TLorentzVector pp = event->GetHMFSParticle(2212)->fP;
// Muon phase space
- //if ((pmu.Vect().Mag() < 250) || cos(pnu.Vect().Angle(pmu.Vect())) < -0.6) {
+ // if ((pmu.Vect().Mag() < 250) || cos(pnu.Vect().Angle(pmu.Vect())) < -0.6) {
// return false;
//}
// Proton phase space
if ((pp.Vect().Mag() < 450) || (cos(pnu.Vect().Angle(pp.Vect())) < 0.4)) {
return false;
}
return true;
}
-}
+} // namespace SignalDef
diff --git a/src/T2K/T2K_SignalDef.h b/src/T2K/T2K_SignalDef.h
index 22325a5..bfc68f5 100644
--- a/src/T2K/T2K_SignalDef.h
+++ b/src/T2K/T2K_SignalDef.h
@@ -1,39 +1,50 @@
// 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/>.
-*******************************************************************************/
+ * 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 T2K_SIGNALDEF_H_SEEN
#define T2K_SIGNALDEF_H_SEEN
#include "SignalDef.h"
namespace SignalDef {
- bool isCC1pip_T2K_H2O(FitEvent *event, double EnuMin, double EnuMax);
- bool isCC1pip_T2K_CH(FitEvent *event, double EnuMin, double EnuMax, bool Michel);
-
- bool isT2K_CC0pi(FitEvent* event, double EnuMin, double EnuMax, bool forwardgoing);
- bool isT2K_CC0piNp(FitEvent* event, double EnuMin, double EnuMax);
- bool isT2K_CC0pi1p(FitEvent* event, double EnuMin, double EnuMax);
- bool isT2K_CC0pi0p(FitEvent* event, double EnuMin, double EnuMax);
- bool isT2K_CC0pi_STV(FitEvent* event, double EnuMin, double EnuMax);
- bool isT2K_CC0pi_1bin(FitEvent* event, double EnuMin, double EnuMax);
- bool isT2K_CC0pi_ifk(FitEvent* event, double EnuMin, double EnuMax);
- bool isT2K_CC0piAnuP0D(FitEvent* event, double EnuMin, double EnuMax);//TN328
-}
+bool isCC1pip_T2K_PRD97_012001(FitEvent *event, double EnuMin, double EnuMax);
+
+enum arxiv1909_03936_PScuts {
+ kMuonFwd = (1 << 0),
+ kMuonHighEff = (1 << 1),
+ kPionFwdHighMom = (1 << 2),
+ kPionHighMom = (1 << 3),
+ kPionHighEff = (1 << 4)
+};
+
+bool isCC1pip_T2K_arxiv1909_03936(FitEvent *event, double EnuMin, double EnuMax,
+ int);
+
+bool isT2K_CC0pi(FitEvent *event, double EnuMin, double EnuMax,
+ bool forwardgoing);
+bool isT2K_CC0piNp(FitEvent *event, double EnuMin, double EnuMax);
+bool isT2K_CC0pi1p(FitEvent *event, double EnuMin, double EnuMax);
+bool isT2K_CC0pi0p(FitEvent *event, double EnuMin, double EnuMax);
+bool isT2K_CC0pi_STV(FitEvent *event, double EnuMin, double EnuMax);
+bool isT2K_CC0pi_1bin(FitEvent *event, double EnuMin, double EnuMax);
+bool isT2K_CC0pi_ifk(FitEvent *event, double EnuMin, double EnuMax);
+bool isT2K_CC0piAnuP0D(FitEvent *event, double EnuMin, double EnuMax); // TN328
+} // namespace SignalDef
#endif