diff --git a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx
index 06646d9..b1b34cb 100644
--- a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx
+++ b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx
@@ -1,353 +1,353 @@
// 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/>.
*******************************************************************************/
// Implementation of 2018 MINERvA numu CC0pi STV analysis
// arxiv 1805.05486.pdf
// Clarence Wret
// cwret@fnal.gov
// Stephen Dolan
// Stephen.Dolan@llr.in2p3.fr
#include "MINERvA_CC0pinp_STV_XSec_1D_nu.h"
#include "MINERvA_SignalDef.h"
//********************************************************************
void MINERvA_CC0pinp_STV_XSec_1D_nu::SetupDataSettings() {
//********************************************************************
// Set Distribution
// See header file for enum and some descriptions
std::string name = fSettings.GetS("name");
if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpmu_nu"))
fDist = kMuonMom;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthmu_nu"))
fDist = kMuonTh;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpprot_nu"))
fDist = kPrMom;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthprot_nu"))
fDist = kPrTh;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpnreco_nu"))
fDist = kNeMom;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddalphat_nu"))
fDist = kDalphaT;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddpt_nu"))
fDist = kDpT;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddphit_nu"))
fDist = kDphiT;
// Location of data, correlation matrices and the titles
std::string titles = "MINERvA_CC0pinp_STV_XSec_1D";
std::string foldername;
std::string distdescript;
// Data release is a single file
std::string rootfile = "MINERvA_1805.05486.root";
fMin = -999;
fMax = 999;
switch (fDist) {
case (kMuonMom):
titles += "pmu";
foldername = "muonmomentum";
distdescript = "Muon momentum in lab frame";
/*
fMin = 2.0;
fMax = 6.0;
*/
fMin = 1.5;
fMax = 10.0;
break;
case (kMuonTh):
titles += "thmu";
foldername = "muontheta";
distdescript = "Muon angle relative neutrino in lab frame";
fMin = 0.0;
fMax = 20.0;
break;
case (kPrMom):
titles += "pprot";
foldername = "protonmomentum";
distdescript = "Proton momentum in lab frame";
// fMin = 0.5;
fMin = 0.45;
fMax = 1.2;
break;
case (kPrTh):
titles += "thprot";
foldername = "protontheta";
distdescript = "Proton angle relative neutrino in lab frame";
fMin = 0.0;
fMax = 70.0;
break;
case (kNeMom):
titles += "pnreco";
foldername = "neutronmomentum";
distdescript = "Neutron momentum in lab frame";
fMin = 0.0;
// fMax = 0.9;
fMax = 2.0;
break;
case (kDalphaT):
foldername = "dalphat";
titles += foldername;
distdescript = "Delta Alpha_T";
fMin = 0.0;
// fMax = 170;
fMax = 180;
break;
case (kDpT):
foldername = "dpt";
titles += foldername;
distdescript = "Delta p_T";
fMin = 0.0;
fMax = 2.0;
break;
case (kDphiT):
foldername = "dphit";
titles += foldername;
distdescript = "Delta phi_T";
fMin = 0.0;
// fMax = 60.0;
fMax = 180.0;
break;
default:
NUIS_ERR(FTL,
"Did not find your specified distribution implemented, exiting");
NUIS_ABORT("You gave " << fName);
}
titles += "_nu";
// All have the same name
std::string dataname = foldername;
// Sample overview ---------------------------------------------------
std::string descrip = distdescript +
"Target: CH \n"
"Flux: MINERvA Forward Horn Current numu ONLY \n"
"Signal: Any event with 1 muon, and 0pi0 in FS, no "
"mesons, at least one proton with: \n"
"1.5GeV < p_mu < 10 GeV\n"
"theta_mu < 20 degrees\n"
"0.45GeV < p_prot < 1.2 GeV\n"
"theta_prot < 70 degrees\n"
"arXiv 1805.05486";
fSettings.SetDescription(descrip);
std::string filename =
GeneralUtils::GetTopLevelDir() +
"/data/MINERvA/CC0pi/CC0pi_STV/MINERvA_1805.05486.root";
// Specify the data
fSettings.SetDataInput(filename + ";" + dataname);
// And the correlations
fSettings.SetCovarInput(filename + ";" + dataname);
// Set titles
fSettings.SetTitle(titles);
};
//********************************************************************
MINERvA_CC0pinp_STV_XSec_1D_nu::MINERvA_CC0pinp_STV_XSec_1D_nu(
nuiskey samplekey) {
//********************************************************************
// A few different distributinos
// Muon momentum, muon angle, proton momentum, proton angle, neutron momentum,
// dalphat, dpt, dphit
// Hard-code the cuts
ProtonMinCut = 450; // MeV
ProtonMaxCut = 1200; // MeV
ProtonThetaCut = 70; // degrees
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
// Load up the data paths and sample descriptions
SetupDataSettings();
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
// No Enu cut
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
// Finalise the settings
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) /
(double(fNEvents) * TotalIntegratedFlux("width"));
// Set the data and covariance matrix
SetDataFromRootFile(fSettings.GetDataInput());
SetCovarianceFromRootFile(fSettings.GetCovarInput());
fSettings.SetXTitle(fDataHist->GetXaxis()->GetTitle());
fSettings.SetYTitle(fDataHist->GetYaxis()->GetTitle());
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
// Data comes in a TList
// Some bins need stripping out because of zero bin content. Why oh why
void MINERvA_CC0pinp_STV_XSec_1D_nu::SetDataFromRootFile(std::string filename) {
//********************************************************************
std::vector<std::string> tempfile = GeneralUtils::ParseToStr(filename, ";");
TFile *File = new TFile(tempfile[0].c_str(), "READ");
// First object is the data
TH1D *temp = (TH1D *)(((TList *)(File->Get(tempfile[1].c_str())))->At(0));
// Garh, some of the data points are zero in the TH1D (WHY?!) so messes with
// the covariance entries to data bins check Skim through the data and check
// for zero bins
std::vector<double> CrossSection;
std::vector<double> Error;
std::vector<double> BinEdges;
int lastbin = 0;
startbin = 0;
for (int i = 0; i < temp->GetXaxis()->GetNbins() + 2; ++i) {
if (temp->GetBinContent(i + 1) > 0 && temp->GetBinLowEdge(i + 1) >= fMin &&
temp->GetBinLowEdge(i + 1) <= fMax) {
if (startbin == 0)
startbin = i;
lastbin = i;
CrossSection.push_back(temp->GetBinContent(i + 1));
BinEdges.push_back(temp->GetXaxis()->GetBinLowEdge(i + 1));
Error.push_back(temp->GetBinError(i + 1));
}
}
BinEdges.push_back(temp->GetXaxis()->GetBinLowEdge(lastbin + 2));
fDataHist = new TH1D((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str(), BinEdges.size() - 1,
&BinEdges[0]);
fDataHist->SetDirectory(0);
for (unsigned int i = 0; i < BinEdges.size() - 1; ++i) {
fDataHist->SetBinContent(i + 1, CrossSection[i]);
fDataHist->SetBinError(i + 1, Error[i]);
}
fDataHist->GetXaxis()->SetTitle(temp->GetXaxis()->GetTitle());
fDataHist->GetYaxis()->SetTitle(temp->GetYaxis()->GetTitle());
fDataHist->SetTitle(temp->GetTitle());
File->Close();
}
//********************************************************************
// Covariance also needs stripping out
// There's padding (two bins...) and overflow (last bin before the two empty
// bins)
void MINERvA_CC0pinp_STV_XSec_1D_nu::SetCovarianceFromRootFile(
std::string filename) {
//********************************************************************
std::vector<std::string> tempfile = GeneralUtils::ParseToStr(filename, ";");
TFile *File = new TFile(tempfile[0].c_str(), "READ");
// First object is the data, second is data with statistical error only, third
// is the covariance matrix
TMatrixDSym *tempcov =
(TMatrixDSym *)((TList *)File->Get(tempfile[1].c_str()))->At(2);
// Count the number of zero entries
int ngood = 0;
int nstart = -1;
int nend = -1;
// Scan through the middle bin and look for entries
int middle = tempcov->GetNrows() / 2;
int nbinsdata = fDataHist->GetXaxis()->GetNbins();
for (int j = 0; j < tempcov->GetNrows(); ++j) {
if ((*tempcov)(middle, j) > 0 && ngood < nbinsdata) {
ngood++;
if (nstart == -1)
nstart = j;
if (j > nend)
nend = j;
}
}
fFullCovar = new TMatrixDSym(ngood);
for (int i = 0; i < fFullCovar->GetNrows(); ++i) {
for (int j = 0; j < fFullCovar->GetNrows(); ++j) {
(*fFullCovar)(i, j) =
(*tempcov)(i + nstart + startbin - 1, j + nstart + startbin - 1);
}
}
(*fFullCovar) *= 1E38 * 1E38;
File->Close();
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
void MINERvA_CC0pinp_STV_XSec_1D_nu::FillEventVariables(FitEvent *event) {
fXVar = -999.99;
// Need a proton and a muon
if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(13) == 0) {
return;
}
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
// Find the highest momentum proton in the event between 450 and 1200 MeV with
// theta_p < 70
TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI));
switch (fDist) {
case (kMuonMom):
fXVar = Pmu.Vect().Mag() / 1000.0;
break;
case (kMuonTh):
fXVar = Pmu.Vect().Angle(Pnu.Vect()) * (180.0 / M_PI);
break;
case (kPrMom):
fXVar = Pprot.Vect().Mag() / 1000.0;
break;
case (kPrTh):
fXVar = Pprot.Vect().Angle(Pnu.Vect()) * (180.0 / M_PI);
break;
// Use Stephen's validated functions
case (kNeMom):
- fXVar = FitUtils::Get_pn_reco_C(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true);
+ fXVar = FitUtils::Get_pn_reco_C_protonps(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true);
break;
case (kDalphaT):
- fXVar = FitUtils::Get_STV_dalphat(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) * (180.0 / M_PI);
+ fXVar = FitUtils::Get_STV_dalphat_protonps(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) * (180.0 / M_PI);
break;
case (kDpT):
- fXVar = FitUtils::Get_STV_dpt(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) / 1000.0;
+ fXVar = FitUtils::Get_STV_dpt_protonps(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) / 1000.0;
break;
case (kDphiT):
- fXVar = FitUtils::Get_STV_dphit(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) * (180.0 / M_PI);
+ fXVar = FitUtils::Get_STV_dphit_protonps(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) * (180.0 / M_PI);
break;
}
return;
};
//********************************************************************
bool MINERvA_CC0pinp_STV_XSec_1D_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isCC0piNp_MINERvA_STV(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
index fbd7423..d322a90 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
@@ -1,94 +1,94 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "T2K_SignalDef.h"
#include "T2K_CC0pinp_STV_XSec_1Ddat_nu.h"
//********************************************************************
T2K_CC0pinp_STV_XSec_1Ddat_nu::T2K_CC0pinp_STV_XSec_1Ddat_nu(
nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0piNp (N>=1) with 450 MeV < p_p < 1 GeV \n"
" p_mu > 250 MeV \n"
" cth_p > 0.4 \n"
" cth_mu > -0.6 \n"
"https://doi.org/10.1103/PhysRevD.98.032003 \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#delta#it{#alpha}_{T} (GeV c^{-1})");
fSettings.SetYTitle(
"#frac{d#sigma}{d#delta#it{#alpha}_{T}} (cm^{2} nucleon^{-1} rads^{-1})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 50.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddat_nu");
// fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() +
// "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.dat");
fSettings.SetDataInput(FitPar::GetDataBase() +
"/T2K/CC0pi/STV/datResults.root;Result");
fSettings.SetCovarInput(FitPar::GetDataBase() +
"/T2K/CC0pi/STV/datResults.root;Correlation_Matrix");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) /
(double(fNEvents) * TotalIntegratedFlux("width"));
// Plot Setup -------------------------------------------------------
// SetDataFromTextFile( fSettings.GetDataInput() );
// ScaleData(1E-38);
// SetCovarFromDiagonal();
SetDataFromRootFile(fSettings.GetDataInput());
SetCorrelationFromRootFile(fSettings.GetCovarInput());
// Proton cuts: 450 MeV to 1 GeV, cos(theta_proton_neutrino) < 0.4
ProtonMinCut = 450;
ProtonMaxCut = 1000;
ProtonCosThetaCut = 0.4;
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void T2K_CC0pinp_STV_XSec_1Ddat_nu::FillEventVariables(FitEvent *event) {
- fXVar = FitUtils::Get_STV_dalphat(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true);
+ fXVar = FitUtils::Get_STV_dalphat_protonps(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true);
return;
};
//********************************************************************
bool T2K_CC0pinp_STV_XSec_1Ddat_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
index 6f7e260..bded9c4 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
@@ -1,93 +1,93 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "T2K_SignalDef.h"
#include "T2K_CC0pinp_STV_XSec_1Ddphit_nu.h"
//********************************************************************
T2K_CC0pinp_STV_XSec_1Ddphit_nu::T2K_CC0pinp_STV_XSec_1Ddphit_nu(
nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0piNp (N>=1) with 450 MeV < p_p < 1 GeV \n"
" p_mu > 250 MeV \n"
" cth_p > 0.4 \n"
" cth_mu > -0.6 \n"
"https://doi.org/10.1103/PhysRevD.98.032003 \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#delta#it{#phi}_{T} (GeV c^{-1})");
fSettings.SetYTitle(
"#frac{d#sigma}{d#delta#it{#phi}_{T}} (cm^{2} nucleon^{-1} rads^{-1})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 50.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddphit_nu");
// fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() +
// "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.dat");
fSettings.SetDataInput(FitPar::GetDataBase() +
"/T2K/CC0pi/STV/dphitResults.root;Result");
fSettings.SetCovarInput(
FitPar::GetDataBase() +
"/T2K/CC0pi/STV/dphitResults.root;Correlation_Matrix");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) /
(double(fNEvents) * TotalIntegratedFlux("width"));
// Plot Setup -------------------------------------------------------
// SetDataFromTextFile( fSettings.GetDataInput() );
// ScaleData(1E-38);
// SetCovarFromDiagonal();
SetDataFromRootFile(fSettings.GetDataInput());
SetCorrelationFromRootFile(fSettings.GetCovarInput());
ProtonMinCut = 450;
ProtonMaxCut = 1000;
ProtonCosThetaCut = 0.4;
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void T2K_CC0pinp_STV_XSec_1Ddphit_nu::FillEventVariables(FitEvent *event) {
- fXVar = FitUtils::Get_STV_dphit(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true);
+ fXVar = FitUtils::Get_STV_dphit_protonps(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true);
return;
};
//********************************************************************
bool T2K_CC0pinp_STV_XSec_1Ddphit_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
index 1f463ab..b44be72 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
@@ -1,93 +1,93 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "T2K_SignalDef.h"
#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h"
//********************************************************************
T2K_CC0pinp_STV_XSec_1Ddpt_nu::T2K_CC0pinp_STV_XSec_1Ddpt_nu(
nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0piNp (N>=1) with 450 MeV < p_p < 1 GeV \n"
" p_mu > 250 MeV \n"
" cth_p > 0.4 \n"
" cth_mu > -0.6 \n"
"https://doi.org/10.1103/PhysRevD.98.032003 \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#delta#it{p}_{T} (GeV c^{-1})");
fSettings.SetYTitle(
"#frac{d#sigma}{d#delta#it{p}_{T}} (cm^{2} nucleon^{-1} GeV^{-1} c)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 50.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddpt_nu");
// fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() +
// "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.dat");
fSettings.SetDataInput(FitPar::GetDataBase() +
"/T2K/CC0pi/STV/dptResults.root;Result");
fSettings.SetCovarInput(FitPar::GetDataBase() +
"/T2K/CC0pi/STV/dptResults.root;Correlation_Matrix");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) /
(double(fNEvents) * TotalIntegratedFlux("width"));
// Plot Setup -------------------------------------------------------
// SetDataFromTextFile( fSettings.GetDataInput() );
// SetCovarFromDiagonal();
// ScaleData(1E-38);
SetDataFromRootFile(fSettings.GetDataInput());
SetCorrelationFromRootFile(fSettings.GetCovarInput());
// SetCovarianceFromRootFile(fSettings.GetCovarInput() );
ProtonMinCut = 450;
ProtonMaxCut = 1000;
ProtonCosThetaCut = 0.4;
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void T2K_CC0pinp_STV_XSec_1Ddpt_nu::FillEventVariables(FitEvent *event) {
- fXVar = FitUtils::Get_STV_dpt(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true) / 1000.0;
+ fXVar = FitUtils::Get_STV_dpt_protonps(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true) / 1000.0;
return;
};
//********************************************************************
bool T2K_CC0pinp_STV_XSec_1Ddpt_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax);
}
diff --git a/src/Utils/FitUtils.cxx b/src/Utils/FitUtils.cxx
index e658bae..5e9d661 100644
--- a/src/Utils/FitUtils.cxx
+++ b/src/Utils/FitUtils.cxx
@@ -1,1098 +1,1158 @@
// 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"
/*
MISC Functions
*/
//********************************************************************
double *FitUtils::GetArrayFromMap(std::vector<std::string> invals,
std::map<std::string, double> inmap) {
//********************************************************************
double *outarr = new double[invals.size()];
int count = 0;
for (size_t i = 0; i < invals.size(); i++) {
outarr[count++] = inmap[invals.at(i)];
}
return outarr;
}
/*
MISC Event
*/
//********************************************************************
// Returns the kinetic energy of a particle in GeV
double FitUtils::T(TLorentzVector part) {
//********************************************************************
double E_part = part.E() / 1000.;
double p_part = part.Vect().Mag() / 1000.;
double m_part = sqrt(E_part * E_part - p_part * p_part);
double KE_part = E_part - m_part;
return KE_part;
};
//********************************************************************
// Returns the momentum of a particle in GeV
double FitUtils::p(TLorentzVector part) {
//********************************************************************
double p_part = part.Vect().Mag() / 1000.;
return p_part;
};
double FitUtils::p(FitParticle *part) { return FitUtils::p(part->fP); };
//********************************************************************
// Returns the angle between two particles in radians
double FitUtils::th(TLorentzVector part1, TLorentzVector part2) {
//********************************************************************
double th = part1.Vect().Angle(part2.Vect());
return th;
};
double FitUtils::th(FitParticle *part1, FitParticle *part2) {
return FitUtils::th(part1->fP, part2->fP);
};
// T2K CC1pi+ helper functions
//
//********************************************************************
// Returns the angle between q3 and the pion defined in Raquel's CC1pi+ on CH
// paper
// Uses "MiniBooNE formula" for Enu, here called EnuCC1pip_T2K_MB
//********************************************************************
double FitUtils::thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
// Want this in GeV
TVector3 p_mu = pmu.Vect() * (1. / 1000.);
// Get the reconstructed Enu
// We are not using Michel e sample, so we have pion kinematic information
double Enu = EnuCC1piprec(pnu, pmu, ppi, true);
// Get neutrino unit direction, multiply by reconstructed Enu
TVector3 p_nu = pnu.Vect() * (1. / (pnu.Vect().Mag())) * Enu;
TVector3 p_pi = ppi.Vect() * (1. / 1000.);
// This is now in GeV
TVector3 q3 = (p_nu - p_mu);
// Want this in GeV
double th_q3_pi = q3.Angle(p_pi);
return th_q3_pi;
}
//********************************************************************
// Returns the q3 defined in Raquel's CC1pi+ on CH paper
// Uses "MiniBooNE formula" for Enu
//********************************************************************
double FitUtils::q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
// Can't use the true Enu here; need to reconstruct it
// We do have Michel e- here so reconstruct Enu by "MiniBooNE formula" without
// pion kinematics
// The bool false refers to that we don't have pion kinematics
// Last bool refers to if we have pion kinematic information or not
double Enu = EnuCC1piprec(pnu, pmu, ppi, false);
TVector3 p_mu = pmu.Vect() * (1. / 1000.);
TVector3 p_nu = pnu.Vect() * (1. / pnu.Vect().Mag()) * Enu;
double q3 = (p_nu - p_mu).Mag();
return q3;
}
//********************************************************************
// Returns the W reconstruction from Raquel CC1pi+ CH thesis
// Uses the MiniBooNE formula Enu
//********************************************************************
double FitUtils::WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
double E_mu = pmu.E() / 1000.;
double p_mu = pmu.Vect().Mag() / 1000.;
double E_nu = EnuCC1piprec(pnu, pmu, ppi, false);
double a1 = (E_nu + PhysConst::mass_neutron) - E_mu;
double a2 = E_nu - p_mu;
double wrec = sqrt(a1 * a1 - a2 * a2);
return wrec;
}
//********************************************************
double FitUtils::ProtonQ2QErec(double pE, double binding) {
//********************************************************
- const double V = binding / 1000.; // binding potential
- const double mn = PhysConst::mass_neutron; // neutron mass
- const double mp = PhysConst::mass_proton; // proton mass
- const double mn_eff = mn - V; // effective proton mass
+ const double V = binding / 1000.; // binding potential
+ const double mn = PhysConst::mass_neutron; // neutron mass
+ const double mp = PhysConst::mass_proton; // proton mass
+ const double mn_eff = mn - V; // effective proton mass
const double pki = (pE / 1000.0) - mp;
double q2qe = mn_eff * mn_eff - mp * mp + 2 * mn_eff * (pki + mp - mn_eff);
return q2qe;
};
//********************************************************************
double FitUtils::EnuQErec(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
- const double V = binding / 1000.; // binding potential
- const double mn = PhysConst::mass_neutron; // neutron mass
- const double mp = PhysConst::mass_proton; // proton mass
+ const double V = binding / 1000.; // binding potential
+ const double mn = PhysConst::mass_neutron; // neutron mass
+ const double mp = PhysConst::mass_proton; // proton mass
double mN_eff = mn - V;
double mN_oth = mp;
if (!neutrino) {
mN_eff = mp - V;
mN_oth = mn;
}
double el = pmu.E() / 1000.;
- double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
- double ml = sqrt(el * el - pl * pl); // lepton mass
+ double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
+ double ml = sqrt(el * el - pl * pl); // lepton mass
double rEnu =
(2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) /
(2 * (mN_eff - el + pl * costh));
return rEnu;
};
// Another good old helper function
-double FitUtils::EnuQErec(TLorentzVector pmu, TLorentzVector pnu, double binding, bool neutrino) {
+double FitUtils::EnuQErec(TLorentzVector pmu, TLorentzVector pnu,
+ double binding, bool neutrino) {
return EnuQErec(pmu, cos(pnu.Vect().Angle(pmu.Vect())), binding, neutrino);
}
-double FitUtils::Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino) {
+double FitUtils::Q2QErec(TLorentzVector pmu, double costh, double binding,
+ bool neutrino) {
double el = pmu.E() / 1000.;
- double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
- double ml = sqrt(el * el - pl * pl); // lepton mass
+ double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
+ double ml = sqrt(el * el - pl * pl); // lepton mass
double rEnu = EnuQErec(pmu, costh, binding, neutrino);
double q2 = -ml * ml + 2. * rEnu * (el - pl * costh);
return q2;
};
-double FitUtils::Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino) {
- double q2qe = Q2QErec(Pmu, cos(Pnu.Vect().Angle(Pmu.Vect())), binding, neutrino);
+double FitUtils::Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding,
+ bool neutrino) {
+ double q2qe =
+ Q2QErec(Pmu, cos(Pnu.Vect().Angle(Pmu.Vect())), binding, neutrino);
return q2qe;
}
double FitUtils::EnuQErec(double pl, double costh, double binding,
bool neutrino) {
- if (pl < 0) return 0.; // Make sure nobody is silly
+ if (pl < 0)
+ return 0.; // Make sure nobody is silly
double mN_eff = PhysConst::mass_neutron - binding / 1000.;
double mN_oth = PhysConst::mass_proton;
if (!neutrino) {
mN_eff = PhysConst::mass_proton - binding / 1000.;
mN_oth = PhysConst::mass_neutron;
}
double ml = PhysConst::mass_muon;
double el = sqrt(pl * pl + ml * ml);
double rEnu =
(2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) /
(2 * (mN_eff - el + pl * costh));
return rEnu;
};
double FitUtils::Q2QErec(double pl, double costh, double binding,
bool neutrino) {
- if (pl < 0) return 0.; // Make sure nobody is silly
+ if (pl < 0)
+ return 0.; // Make sure nobody is silly
double ml = PhysConst::mass_muon;
double el = sqrt(pl * pl + ml * ml);
double rEnu = EnuQErec(pl, costh, binding, neutrino);
double q2 = -ml * ml + 2. * rEnu * (el - pl * costh);
return q2;
};
//********************************************************************
// Reconstructs Enu for CC1pi0
// Very similar for CC1pi+ reconstruction
double FitUtils::EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0) {
//********************************************************************
double E_mu = pmu.E() / 1000;
double p_mu = pmu.Vect().Mag() / 1000;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double E_pi0 = ppi0.E() / 1000;
double p_pi0 = ppi0.Vect().Mag() / 1000;
double m_pi0 = sqrt(E_pi0 * E_pi0 - p_pi0 * p_pi0);
double th_nu_pi0 = pnu.Vect().Angle(ppi0.Vect());
- const double m_n = PhysConst::mass_neutron; // neutron mass
- const double m_p = PhysConst::mass_proton; // proton mass
+ const double m_n = PhysConst::mass_neutron; // neutron mass
+ const double m_p = PhysConst::mass_proton; // proton mass
double th_pi0_mu = ppi0.Vect().Angle(pmu.Vect());
double rEnu = (m_mu * m_mu + m_pi0 * m_pi0 + m_n * m_n - m_p * m_p -
2 * m_n * (E_pi0 + E_mu) + 2 * E_pi0 * E_mu -
2 * p_pi0 * p_mu * cos(th_pi0_mu)) /
(2 * (E_pi0 + E_mu - p_pi0 * cos(th_nu_pi0) -
p_mu * cos(th_nu_mu) - m_n));
return rEnu;
};
//********************************************************************
// Reconstruct Q2 for CC1pi0
// Beware: uses true Enu, not reconstructed Enu
double FitUtils::Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0) {
//********************************************************************
- double E_mu = pmu.E() / 1000.; // energy of lepton in GeV
- double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton
- double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass
+ double E_mu = pmu.E() / 1000.; // energy of lepton in GeV
+ double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton
+ double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// double rEnu = EnuCC1pi0rec(pnu, pmu, ppi0); //reconstructed neutrino energy
// Use true neutrino energy
double rEnu = pnu.E() / 1000.;
double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu));
return q2;
};
//********************************************************************
// Reconstruct Enu for CC1pi+
// pionInfo reflects if we're using pion kinematics or not
// In T2K CC1pi+ CH the Michel tag is used for pion in which pion kinematic info
// is lost and Enu is reconstructed without pion kinematics
double FitUtils::EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi, bool pionInfo) {
//********************************************************************
double E_mu = pmu.E() / 1000.;
double p_mu = pmu.Vect().Mag() / 1000.;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double E_pi = ppi.E() / 1000.;
double p_pi = ppi.Vect().Mag() / 1000.;
double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi);
- const double m_n = PhysConst::mass_neutron; // neutron/proton mass
+ const double m_n = PhysConst::mass_neutron; // neutron/proton mass
// should really take proton mass for proton interaction, neutron for neutron
// interaction. However, difference is pretty much negligable here!
// need this angle too
double th_nu_pi = pnu.Vect().Angle(ppi.Vect());
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double th_pi_mu = ppi.Vect().Angle(pmu.Vect());
double rEnu = -999;
// If experiment doesn't have pion kinematic information (T2K CC1pi+ CH
// example of this; Michel e sample has no directional information on pion)
// We'll need to modify the reconstruction variables
if (!pionInfo) {
// Assumes that pion angle contribution contributes net zero
// Also assumes the momentum of reconstructed pions when using Michel
// electrons is 130 MeV
// So need to redefine E_pi and p_pi
// It's a little unclear what happens to the pmu . ppi term (4-vector); do
// we include the angular contribution?
// This below doesn't
th_nu_pi = M_PI / 2.;
p_pi = 0.130;
E_pi = sqrt(p_pi * p_pi + m_pi * m_pi);
// rEnu = (m_mu*m_mu + m_pi*m_pi - 2*m_n*(E_mu + E_pi) + 2*E_mu*E_pi -
// 2*p_mu*p_pi) / (2*(E_mu + E_pi - p_mu*cos(th_nu_mu) - m_n));
}
rEnu =
(m_mu * m_mu + m_pi * m_pi - 2 * m_n * (E_pi + E_mu) + 2 * E_pi * E_mu -
2 * p_pi * p_mu * cos(th_pi_mu)) /
(2 * (E_pi + E_mu - p_pi * cos(th_nu_pi) - p_mu * cos(th_nu_mu) - m_n));
return rEnu;
};
//********************************************************************
// Reconstruct neutrino energy from outgoing particles; will differ from the
// actual neutrino energy. Here we use assumption of a Delta resonance
double FitUtils::EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu) {
//********************************************************************
const double m_Delta =
- PhysConst::mass_delta; // PDG value for Delta mass in GeV
- const double m_n = PhysConst::mass_neutron; // neutron/proton mass
+ PhysConst::mass_delta; // PDG value for Delta mass in GeV
+ const double m_n = PhysConst::mass_neutron; // neutron/proton mass
// should really take proton mass for proton interaction, neutron for neutron
// interaction. However, difference is pretty much negligable here!
double E_mu = pmu.E() / 1000;
double p_mu = pmu.Vect().Mag() / 1000;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double rEnu = (m_Delta * m_Delta - m_n * m_n - m_mu * m_mu + 2 * m_n * E_mu) /
(2 * (m_n - E_mu + p_mu * cos(th_nu_mu)));
return rEnu;
};
// MOVE TO T2K UTILS!
//********************************************************************
// Reconstruct Enu using "extended MiniBooNE" as defined in Raquel's T2K TN
//
// Supposedly includes pion direction and binding energy of target nucleon
// I'm not convinced (yet), maybe
double FitUtils::EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
//********************************************************************
// Unit vector for neutrino momentum
TVector3 p_nu_vect_unit = pnu.Vect() * (1. / pnu.E());
double E_mu = pmu.E() / 1000.;
TVector3 p_mu_vect = pmu.Vect() * (1. / 1000.);
double E_pi = ppi.E() / 1000.;
TVector3 p_pi_vect = ppi.Vect() * (1. / 1000.);
- double E_bind = 25. / 1000.;
+ double E_bind = 25. / 1000.;
double m_p = PhysConst::mass_proton;
// Makes life a little easier, gonna square this one
double a1 = m_p - E_bind - E_mu - E_pi;
// Just gets the magnitude square of the muon and pion momentum vectors
double a2 = (p_mu_vect + p_pi_vect).Mag2();
// Gets the somewhat complicated scalar product between neutrino and
// (p_mu+p_pi), scaled to Enu
double a3 = p_nu_vect_unit * (p_mu_vect + p_pi_vect);
double rEnu =
(m_p * m_p - a1 * a1 + a2) / (2 * (m_p - E_bind - E_mu - E_pi + a3));
return rEnu;
}
//********************************************************************
// Reconstructed Q2 for CC1pi+
//
// enuType describes how the neutrino energy is reconstructed
// 0 uses true neutrino energy; case for MINERvA and MiniBooNE
// 1 uses "extended MiniBooNE" reconstructed (T2K CH)
// 2 uses "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = true) (T2K CH)
// "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = false, the
// case for T2K when using Michel tag) (T2K CH)
// 3 uses Delta for reconstruction (T2K CH)
double FitUtils::Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi, int enuType, bool pionInfo) {
//********************************************************************
- double E_mu = pmu.E() / 1000.; // energy of lepton in GeV
- double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton
- double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass
+ double E_mu = pmu.E() / 1000.; // energy of lepton in GeV
+ double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton
+ double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// Different ways of reconstructing the neutrino energy; default is just to
// use the truth (case 0)
double rEnu = -999;
switch (enuType) {
- case 0: // True neutrino energy, default; this is the case for when Q2
- // defined as Q2 true (MiniBooNE, MINERvA)
- rEnu = pnu.E() / 1000.;
- break;
- case 1: // Extended MiniBooNE reconstructed, as defined by Raquel's CC1pi+
- // CH T2K analysis
- rEnu = EnuCC1piprec_T2K_eMB(pnu, pmu, ppi);
- break;
- case 2: // MiniBooNE reconstructed, as defined by MiniBooNE and Raquel's
- // CC1pi+ CH T2K analysis and Linda's CC1pi+ H2O
- // Can have this with and without pion info, depending on if Michel tag
- // used (Raquel)
- rEnu = EnuCC1piprec(pnu, pmu, ppi, pionInfo);
- break;
- case 3:
- rEnu = EnuCC1piprecDelta(pnu, pmu);
- break;
-
- } // No need for default here since default value of enuType = 0, defined in
- // header
+ case 0: // True neutrino energy, default; this is the case for when Q2
+ // defined as Q2 true (MiniBooNE, MINERvA)
+ rEnu = pnu.E() / 1000.;
+ break;
+ case 1: // Extended MiniBooNE reconstructed, as defined by Raquel's CC1pi+
+ // CH T2K analysis
+ rEnu = EnuCC1piprec_T2K_eMB(pnu, pmu, ppi);
+ break;
+ case 2: // MiniBooNE reconstructed, as defined by MiniBooNE and Raquel's
+ // CC1pi+ CH T2K analysis and Linda's CC1pi+ H2O
+ // Can have this with and without pion info, depending on if Michel tag
+ // used (Raquel)
+ rEnu = EnuCC1piprec(pnu, pmu, ppi, pionInfo);
+ break;
+ case 3:
+ rEnu = EnuCC1piprecDelta(pnu, pmu);
+ break;
+
+ } // No need for default here since default value of enuType = 0, defined in
+ // header
double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu));
return q2;
};
//********************************************************************
// Returns the reconstructed W from a nucleon and an outgoing pion
//
// Could do this a lot more clever (pp + ppi).Mag() would do the job, but this
// would be less instructive
//********************************************************************
double FitUtils::MpPi(TLorentzVector pp, TLorentzVector ppi) {
double E_p = pp.E();
double p_p = pp.Vect().Mag();
double m_p = sqrt(E_p * E_p - p_p * p_p);
double E_pi = ppi.E();
double p_pi = ppi.Vect().Mag();
double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi);
double th_p_pi = pp.Vect().Angle(ppi.Vect());
// fairly easy thing to derive since bubble chambers measure the proton!
double invMass = sqrt(m_p * m_p + m_pi * m_pi + 2 * E_p * E_pi -
2 * p_pi * p_p * cos(th_p_pi));
return invMass;
};
//********************************************************
// Reconstruct the hadronic mass using neutrino and muon
// Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be
// reconstructed Enu; so gives reconstructed Wrec which most of the time isn't
// used!
//
// Only MINERvA uses this so far; and the Enu is Enu_true
// If we want W_true need to take initial state nucleon motion into account
// Return value is in MeV!!!
double FitUtils::Wrec(TLorentzVector pnu, TLorentzVector pmu) {
//********************************************************
double E_mu = pmu.E();
double p_mu = pmu.Vect().Mag();
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// The factor of 1000 is necessary for downstream functions
const double m_p = PhysConst::mass_proton * 1000;
// MINERvA cut on W_exp which is tuned to W_true; so use true Enu from
// generators
double E_nu = pnu.E();
double w_rec = sqrt(m_p * m_p + m_mu * m_mu - 2 * m_p * E_mu +
2 * E_nu * (m_p - E_mu + p_mu * cos(th_nu_mu)));
return w_rec;
};
//********************************************************
// Reconstruct the true hadronic mass using the initial state and muon
// Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be
// reconstructed Enu; so gives reconstructed Wrec which most of the time isn't
// used!
//
// No one seems to use this because it's fairly MC dependent!
// Return value is in MeV!!!
double FitUtils::Wtrue(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector pnuc) {
//********************************************************
// Could simply do the TLorentzVector operators here but this is more
// instructive?
// ... and prone to errors ...
double E_mu = pmu.E();
double p_mu = pmu.Vect().Mag();
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double E_nuc = pnuc.E();
double p_nuc = pnuc.Vect().Mag();
double m_nuc = sqrt(E_nuc * E_nuc - p_nuc * p_nuc);
double th_nuc_mu = pmu.Vect().Angle(pnuc.Vect());
double th_nu_nuc = pnu.Vect().Angle(pnuc.Vect());
double E_nu = pnu.E();
double w_rec = sqrt(m_nuc * m_nuc + m_mu * m_mu - 2 * E_nu * E_mu +
2 * E_nu * p_mu * cos(th_nu_mu) - 2 * E_nuc * E_mu +
2 * p_nuc * p_mu * cos(th_nuc_mu) + 2 * E_nu * E_nuc -
2 * E_nu * p_nuc * cos(th_nu_nuc));
return w_rec;
};
double FitUtils::SumKE_PartVect(std::vector<FitParticle *> const fps) {
double sum = 0.0;
for (size_t p_it = 0; p_it < fps.size(); ++p_it) {
sum += fps[p_it]->KE();
}
return sum;
}
double FitUtils::SumTE_PartVect(std::vector<FitParticle *> const fps) {
double sum = 0.0;
for (size_t p_it = 0; p_it < fps.size(); ++p_it) {
sum += fps[p_it]->E();
}
return sum;
}
/*
E Recoil
*/
double FitUtils::GetErecoil_TRUE(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
- if (!event->PartInfo(i)->fIsAlive) continue;
- if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
+ if (!event->PartInfo(i)->fIsAlive)
+ continue;
+ if (event->PartInfo(i)->fNEUTStatusCode != 0)
+ continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1)
continue;
// Add Up KE of protons and TE of everything else
if (event->PartInfo(i)->fPID == 2212 || event->PartInfo(i)->fPID == 2112) {
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
} else {
Erecoil += event->PartInfo(i)->fP.E();
}
}
return Erecoil;
}
double FitUtils::GetErecoil_CHARGED(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
- if (!event->PartInfo(i)->fIsAlive) continue;
- if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
+ if (!event->PartInfo(i)->fIsAlive)
+ continue;
+ if (event->PartInfo(i)->fNEUTStatusCode != 0)
+ continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1)
continue;
// Skip Neutral particles
if (event->PartInfo(i)->fPID == 2112 || event->PartInfo(i)->fPID == 111 ||
event->PartInfo(i)->fPID == 22)
continue;
// Add Up KE of protons and TE of everything else
if (event->PartInfo(i)->fPID == 2212) {
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
} else {
Erecoil += event->PartInfo(i)->fP.E();
}
}
return Erecoil;
}
// MOVE TO MINERVA Utils!
double FitUtils::GetErecoil_MINERvA_LowRecoil(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
- if (!event->PartInfo(i)->fIsAlive) continue;
- if (event->PartInfo(i)->Status() != kFinalState) continue;
+ if (!event->PartInfo(i)->fIsAlive)
+ continue;
+ if (event->PartInfo(i)->Status() != kFinalState)
+ continue;
// Skip Lepton
- if (abs(event->PartInfo(i)->fPID) == 13) continue;
+ if (abs(event->PartInfo(i)->fPID) == 13)
+ continue;
// Skip Neutrons particles
- if (event->PartInfo(i)->fPID == 2112) continue;
+ if (event->PartInfo(i)->fPID == 2112)
+ continue;
int PID = event->PartInfo(i)->fPID;
// KE of Protons and charged pions
if (PID == 2212 or PID == 211 or PID == -211) {
// Erecoil += FitUtils::T(event->PartInfo(i)->fP);
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
// Total Energy of non-neutrons
// } else if (PID != 2112 and PID < 999 and PID != 22 and abs(PID) !=
// 14) {
} else if (PID == 111 || PID == 11 || PID == -11 || PID == 22) {
Erecoil += (event->PartInfo(i)->fP.E());
}
}
return Erecoil;
}
// MOVE TO MINERVA Utils!
-// The alternative Eavailble definition takes true q0 and subtracts the kinetic energy of neutrons and pion masses
-// returns in MeV
+// The alternative Eavailble definition takes true q0 and subtracts the kinetic
+// energy of neutrons and pion masses returns in MeV
double FitUtils::Eavailable(FitEvent *event) {
double Eav = 0.0;
// Now take q0 and subtract Eav
double q0 = event->GetBeamPart()->fP.E();
// Get the pdg of incoming neutrino
int ISPDG = event->GetBeamPartPDG();
// For CC
- if (event->IsCC()) q0 -= event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.E();
- else q0 -= event->GetHMFSParticle(ISPDG)->fP.E();
+ if (event->IsCC())
+ q0 -= event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.E();
+ else
+ q0 -= event->GetHMFSParticle(ISPDG)->fP.E();
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
- if (!event->PartInfo(i)->fIsAlive) continue;
- if (event->PartInfo(i)->Status() != kFinalState) continue;
+ if (!event->PartInfo(i)->fIsAlive)
+ continue;
+ if (event->PartInfo(i)->Status() != kFinalState)
+ continue;
int PID = event->PartInfo(i)->fPID;
// Neutrons
if (PID == 2112) {
// Adding kinetic energy of neutron
- Eav += FitUtils::T(event->PartInfo(i)->fP)*1000.;
+ Eav += FitUtils::T(event->PartInfo(i)->fP) * 1000.;
// All pion masses
} else if (abs(PID) == 211 || PID == 111) {
Eav += event->PartInfo(i)->fP.M();
}
}
- return q0-Eav;
+ return q0 - Eav;
}
TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) {
TVector3 pnUnit = planarNormal.Unit();
double inpProjectPN = inp.Dot(pnUnit);
TVector3 InPlanarInput = inp - (inpProjectPN * pnUnit);
return InPlanarInput;
}
TVector3 GetUnitVectorInTPlane(const TVector3 &inp,
const TVector3 &planarNormal) {
return GetVectorInTPlane(inp, planarNormal).Unit();
}
Double_t GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal, bool PiMinus = false) {
TVector3 V_lepton_T = GetUnitVectorInTPlane(V_lepton, Normal);
TVector3 V_other_T = GetUnitVectorInTPlane(V_other, Normal);
return PiMinus ? acos(V_lepton_T.Dot(V_other_T))
: (M_PI - acos(V_lepton_T.Dot(V_other_T)));
}
TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal) {
TVector3 V_lepton_T = GetVectorInTPlane(V_lepton, Normal);
TVector3 V_other_T = GetVectorInTPlane(V_other, Normal);
return V_lepton_T + V_other_T;
}
Double_t GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal, bool PiMinus = false) {
TVector3 DeltaPT = GetDeltaPT(V_lepton, V_other, Normal);
return GetDeltaPhiT(V_lepton, DeltaPT, Normal, PiMinus);
}
-double FitUtils::Get_STV_dpt(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut, int ISPDG, bool Is0pi) {
+double FitUtils::Get_STV_dpt_protonps(FitEvent *event, double ProtonMinCut,
+ double ProtonMaxCut,
+ double ProtonCosThetaCut, int ISPDG,
+ bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
- // Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
- TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
+ // Find the highest momentum proton in the event between ProtonMinCut and
+ // ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = FitUtils::GetProtonInRange(
+ event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
+
// Get highest momentum proton in allowed proton range
TVector3 HadronP = Pprot.Vect();
// If we don't have a CC0pi signal definition we also add in pion momentum
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
// Count up pion momentum
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaPT(LeptonP, HadronP, NuP).Mag();
}
-double FitUtils::Get_STV_dphit(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut, int ISPDG, bool Is0pi) {
+double FitUtils::Get_STV_dphit_protonps(FitEvent *event, double ProtonMinCut,
+ double ProtonMaxCut,
+ double ProtonCosThetaCut, int ISPDG,
+ bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
- // Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
- TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
+ // Find the highest momentum proton in the event between ProtonMinCut and
+ // ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = FitUtils::GetProtonInRange(
+ event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
TVector3 HadronP = Pprot.Vect();
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaPhiT(LeptonP, HadronP, NuP);
}
-double FitUtils::Get_STV_dalphat(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut, int ISPDG, bool Is0pi) {
+double FitUtils::Get_STV_dalphat_protonps(FitEvent *event, double ProtonMinCut,
+ double ProtonMaxCut,
+ double ProtonCosThetaCut, int ISPDG,
+ bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
- // Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
- TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
+ // Find the highest momentum proton in the event between ProtonMinCut and
+ // ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = FitUtils::GetProtonInRange(
+ event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
TVector3 HadronP = Pprot.Vect();
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaAlphaT(LeptonP, HadronP, NuP);
}
// As defined in PhysRevC.95.065501
-// Using prescription from arXiv 1805.05486
+// Using prescription from arXiv 1805.05486
// Returns in GeV
-double FitUtils::Get_pn_reco_C(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut, int ISPDG, bool Is0pi) {
+double FitUtils::Get_pn_reco_C_protonps(FitEvent *event, double ProtonMinCut,
+ double ProtonMaxCut,
+ double ProtonCosThetaCut, int ISPDG,
+ bool Is0pi) {
- const double mn = PhysConst::mass_neutron; // neutron mass
- const double mp = PhysConst::mass_proton; // proton mass
+ const double mn = PhysConst::mass_neutron; // neutron mass
+ const double mp = PhysConst::mass_proton; // proton mass
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
- // Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) < ProtonCosThetaCut
- TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
+ // Find the highest momentum proton in the event between ProtonMinCut and
+ // ProtonMaxCut MeV with cos(theta_p) < ProtonCosThetaCut
+ TLorentzVector Pprot = FitUtils::GetProtonInRange(
+ event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
TVector3 HadronP = Pprot.Vect();
- double const el = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E()/1000.;
- double const eh = Pprot.E()/1000.;
+ double const el =
+ event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E() / 1000.;
+ double const eh = Pprot.E() / 1000.;
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP);
- double dptMag = dpt.Mag()/1000.;
+ double dptMag = dpt.Mag() / 1000.;
- double ma = 6*mn + 6*mp - 0.09216; // target mass (E is from PhysRevC.95.065501)
+ double ma =
+ 6 * mn + 6 * mp - 0.09216; // target mass (E is from PhysRevC.95.065501)
double map = ma - mn + 0.02713; // remnant mass
- double pmul = LeptonP.Dot(NuP.Unit())/1000.;
- double phl = HadronP.Dot(NuP.Unit())/1000.;
+ double pmul = LeptonP.Dot(NuP.Unit()) / 1000.;
+ double phl = HadronP.Dot(NuP.Unit()) / 1000.;
- //double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.;
- //double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.;
+ // double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.;
+ // double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.;
double R = ma + pmul + phl - el - eh;
- double dpl = 0.5*R - (map*map + dptMag*dptMag)/(2*R);
- //double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in PhysRevC.95.065501 - gives same result
-
- double pn_reco = sqrt((dptMag*dptMag) + (dpl*dpl));
-
- //std::cout << "Diagnostics: " << std::endl;
- //std::cout << "mn: " << mn << std::endl;
- //std::cout << "ma: " << ma << std::endl;
- //std::cout << "map: " << map << std::endl;
- //std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl;
- //std::cout << "ph: " << HadronP.Mag()/1000. << std::endl;
- //std::cout << "pmul: " << pmul << std::endl;
- //std::cout << "phl: " << phl << std::endl;
- //std::cout << "el: " << el << std::endl;
- //std::cout << "eh: " << eh << std::endl;
- //std::cout << "R: " << R << std::endl;
- //std::cout << "dptMag: " << dptMag << std::endl;
- //std::cout << "dpl: " << dpl << std::endl;
- //std::cout << "pn_reco: " << pn_reco << std::endl;
+ double dpl = 0.5 * R - (map * map + dptMag * dptMag) / (2 * R);
+ // double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in
+ // PhysRevC.95.065501 - gives same result
+
+ double pn_reco = sqrt((dptMag * dptMag) + (dpl * dpl));
+
+ // std::cout << "Diagnostics: " << std::endl;
+ // std::cout << "mn: " << mn << std::endl;
+ // std::cout << "ma: " << ma << std::endl;
+ // std::cout << "map: " << map << std::endl;
+ // std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl;
+ // std::cout << "ph: " << HadronP.Mag()/1000. << std::endl;
+ // std::cout << "pmul: " << pmul << std::endl;
+ // std::cout << "phl: " << phl << std::endl;
+ // std::cout << "el: " << el << std::endl;
+ // std::cout << "eh: " << eh << std::endl;
+ // std::cout << "R: " << R << std::endl;
+ // std::cout << "dptMag: " << dptMag << std::endl;
+ // std::cout << "dpl: " << dpl << std::endl;
+ // std::cout << "pn_reco: " << pn_reco << std::endl;
return pn_reco;
}
-double FitUtils::Get_pn_reco_Ar(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut, int ISPDG, bool Is0pi) {
+double FitUtils::Get_pn_reco_Ar_protonps(FitEvent *event, double ProtonMinCut,
+ double ProtonMaxCut,
+ double ProtonCosThetaCut, int ISPDG,
+ bool Is0pi) {
- const double mn = PhysConst::mass_neutron; // neutron mass
- const double mp = PhysConst::mass_proton; // proton mass
+ const double mn = PhysConst::mass_neutron; // neutron mass
+ const double mp = PhysConst::mass_proton; // proton mass
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
- // Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
- TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
+ // Find the highest momentum proton in the event between ProtonMinCut and
+ // ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = FitUtils::GetProtonInRange(
+ event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
TVector3 HadronP = Pprot.Vect();
- double const el = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E()/1000.;
- double const eh = Pprot.E()/1000.;
+ double const el =
+ event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E() / 1000.;
+ double const eh = Pprot.E() / 1000.;
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP);
- double dptMag = dpt.Mag()/1000.;
+ double dptMag = dpt.Mag() / 1000.;
- //double ma = 6*mn + 6*mp - 0.09216; // target mass (E is from PhysRevC.95.065501)
- //double map = ma - mn + 0.02713; // remnant mass
- double ma = 6*mn + 6*mp - 0.34381; // target mass (E is from PhysRevC.95.065501)
+ // double ma = 6*mn + 6*mp - 0.09216; // target mass (E is from
+ // PhysRevC.95.065501) double map = ma - mn + 0.02713; // remnant mass
+ double ma =
+ 6 * mn + 6 * mp - 0.34381; // target mass (E is from PhysRevC.95.065501)
double map = ma - mn + 0.02713; // remnant mass
- double pmul = LeptonP.Dot(NuP.Unit())/1000.;
- double phl = HadronP.Dot(NuP.Unit())/1000.;
+ double pmul = LeptonP.Dot(NuP.Unit()) / 1000.;
+ double phl = HadronP.Dot(NuP.Unit()) / 1000.;
- //double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.;
- //double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.;
+ // double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.;
+ // double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.;
double R = ma + pmul + phl - el - eh;
- double dpl = 0.5*R - (map*map + dptMag*dptMag)/(2*R);
- //double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in PhysRevC.95.065501 - gives same result
-
- double pn_reco = sqrt((dptMag*dptMag) + (dpl*dpl));
-
- //std::cout << "Diagnostics: " << std::endl;
- //std::cout << "mn: " << mn << std::endl;
- //std::cout << "ma: " << ma << std::endl;
- //std::cout << "map: " << map << std::endl;
- //std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl;
- //std::cout << "ph: " << HadronP.Mag()/1000. << std::endl;
- //std::cout << "pmul: " << pmul << std::endl;
- //std::cout << "phl: " << phl << std::endl;
- //std::cout << "el: " << el << std::endl;
- //std::cout << "eh: " << eh << std::endl;
- //std::cout << "R: " << R << std::endl;
- //std::cout << "dptMag: " << dptMag << std::endl;
- //std::cout << "dpl: " << dpl << std::endl;
- //std::cout << "pn_reco: " << pn_reco << std::endl;
+ double dpl = 0.5 * R - (map * map + dptMag * dptMag) / (2 * R);
+ // double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in
+ // PhysRevC.95.065501 - gives same result
+
+ double pn_reco = sqrt((dptMag * dptMag) + (dpl * dpl));
+
+ // std::cout << "Diagnostics: " << std::endl;
+ // std::cout << "mn: " << mn << std::endl;
+ // std::cout << "ma: " << ma << std::endl;
+ // std::cout << "map: " << map << std::endl;
+ // std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl;
+ // std::cout << "ph: " << HadronP.Mag()/1000. << std::endl;
+ // std::cout << "pmul: " << pmul << std::endl;
+ // std::cout << "phl: " << phl << std::endl;
+ // std::cout << "el: " << el << std::endl;
+ // std::cout << "eh: " << eh << std::endl;
+ // std::cout << "R: " << R << std::endl;
+ // std::cout << "dptMag: " << dptMag << std::endl;
+ // std::cout << "dpl: " << dpl << std::endl;
+ // std::cout << "pn_reco: " << pn_reco << std::endl;
return pn_reco;
}
-// Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
-TLorentzVector FitUtils::GetProtonInRange(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut) {
+// Find the highest momentum proton in the event between ProtonMinCut and
+// ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
+TLorentzVector FitUtils::GetProtonInRange(FitEvent *event, double ProtonMinCut,
+ double ProtonMaxCut,
+ double ProtonCosThetaCut) {
// Get the neutrino
- TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
+ TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
int HMFSProton = 0;
double HighestMomentum = 0.0;
// Get the stack of protons
- std::vector<FitParticle*> Protons = event->GetAllFSProton();
+ std::vector<FitParticle *> Protons = event->GetAllFSProton();
for (size_t i = 0; i < Protons.size(); ++i) {
- if (Protons[i]->p() > ProtonMinCut &&
- Protons[i]->p() < ProtonMaxCut &&
+ if (Protons[i]->p() > ProtonMinCut && Protons[i]->p() < ProtonMaxCut &&
cos(Protons[i]->P3().Angle(Pnu.Vect())) > ProtonCosThetaCut &&
Protons[i]->p() > HighestMomentum) {
HighestMomentum = Protons[i]->p();
HMFSProton = i;
}
}
// Now get the proton
TLorentzVector Pprot = Protons[HMFSProton]->fP;
return Pprot;
}
// Get Cos theta with Adler angles
-double FitUtils::CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) {
+double FitUtils::CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu,
+ TLorentzVector Ppi, TLorentzVector Pprot) {
// Get the "resonance" lorentz vector (pion proton system)
TLorentzVector Pres = Pprot + Ppi;
- // Boost the particles into the resonance rest frame so we can define the x,y,z axis
+ // Boost the particles into the resonance rest frame so we can define the
+ // x,y,z axis
Pnu.Boost(-Pres.BoostVector());
Pmu.Boost(-Pres.BoostVector());
Ppi.Boost(-Pres.BoostVector());
// The z-axis is defined as the axis of three-momentum transfer, \vec{k}
// Also unit normalise the axis
- TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag()));
+ TVector3 zAxis =
+ (Pnu.Vect() - Pmu.Vect()) * (1.0 / ((Pnu.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(Ppi.Vect().Angle(zAxis));
return costhAdler;
}
// Get phi with Adler angles, a bit more complicated...
-double FitUtils::PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) {
+double FitUtils::PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu,
+ TLorentzVector Ppi, TLorentzVector Pprot) {
// Get the "resonance" lorentz vector (pion proton system)
TLorentzVector Pres = Pprot + Ppi;
- // Boost the particles into the resonance rest frame so we can define the x,y,z axis
+ // Boost the particles into the resonance rest frame so we can define the
+ // x,y,z axis
Pnu.Boost(-Pres.BoostVector());
Pmu.Boost(-Pres.BoostVector());
Ppi.Boost(-Pres.BoostVector());
// The z-axis is defined as the axis of three-momentum transfer, \vec{k}
// Also unit normalise the axis
- TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag()));
+ TVector3 zAxis =
+ (Pnu.Vect() - Pmu.Vect()) * (1.0 / ((Pnu.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 = Pnu.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());
// Project the pion on to x and y axes
double x = Ppi.Vect().Dot(xAxis);
double y = Ppi.Vect().Dot(yAxis);
- double newphi = atan2(y, x)*(180./M_PI);
+ double newphi = atan2(y, x) * (180. / M_PI);
// Convert negative angles to positive
- if (newphi < 0.0) newphi += 360.0;
+ if (newphi < 0.0)
+ newphi += 360.0;
return newphi;
}
-
//********************************************************************
double FitUtils::ppInfK(TLorentzVector pmu, double costh, double binding,
- bool neutrino) {
+ bool neutrino) {
//********************************************************************
// Convert all values to GeV
- //const double V = binding / 1000.; // binding potential
- //const double mn = PhysConst::mass_neutron; // neutron mass
- const double mp = PhysConst::mass_proton; // proton mass
+ // const double V = binding / 1000.; // binding potential
+ // const double mn = PhysConst::mass_neutron; // neutron mass
+ const double mp = PhysConst::mass_proton; // proton mass
double el = pmu.E() / 1000.;
- //double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
+ // double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double enu = EnuQErec(pmu, costh, binding, neutrino);
double ep_inf = enu - el + mp;
double pp_inf = sqrt(ep_inf * ep_inf - mp * mp);
return pp_inf;
};
//********************************************************************
TVector3 FitUtils::tppInfK(TLorentzVector pmu, double costh, double binding,
- bool neutrino) {
+ bool neutrino) {
//********************************************************************
// Convert all values to GeV
- //const double V = binding / 1000.; // binding potential
- //const double mn = PhysConst::mass_neutron; // neutron mass
- //const double mp = PhysConst::mass_proton; // proton mass
+ // const double V = binding / 1000.; // binding potential
+ // const double mn = PhysConst::mass_neutron; // neutron mass
+ // const double mp = PhysConst::mass_proton; // proton mass
double pl_x = pmu.X() / 1000.;
double pl_y = pmu.Y() / 1000.;
- double pl_z= pmu.Z() / 1000.;
+ double pl_z = pmu.Z() / 1000.;
double enu = EnuQErec(pmu, costh, binding, neutrino);
- TVector3 tpp_inf(-pl_x, -pl_y, -pl_z+enu);
+ TVector3 tpp_inf(-pl_x, -pl_y, -pl_z + enu);
return tpp_inf;
};
//********************************************************************
double FitUtils::cthpInfK(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
- //const double V = binding / 1000.; // binding potential
- //const double mn = PhysConst::mass_neutron; // neutron mass
- const double mp = PhysConst::mass_proton; // proton mass
+ // const double V = binding / 1000.; // binding potential
+ // const double mn = PhysConst::mass_neutron; // neutron mass
+ const double mp = PhysConst::mass_proton; // proton mass
double el = pmu.E() / 1000.;
- double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
+ double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double enu = EnuQErec(pmu, costh, binding, neutrino);
double ep_inf = enu - el + mp;
double pp_inf = sqrt(ep_inf * ep_inf - mp * mp);
- double cth_inf = (enu*enu + pp_inf*pp_inf - pl*pl)/(2*enu*pp_inf);
+ double cth_inf = (enu * enu + pp_inf * pp_inf - pl * pl) / (2 * enu * pp_inf);
return cth_inf;
};
diff --git a/src/Utils/FitUtils.h b/src/Utils/FitUtils.h
index cc9326c..1ba9e85 100644
--- a/src/Utils/FitUtils.h
+++ b/src/Utils/FitUtils.h
@@ -1,193 +1,256 @@
// 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 FITUTILS_H_SEEN
#define FITUTILS_H_SEEN
-#include <math.h>
-#include <stdlib.h>
-#include <unistd.h>
#include <ctime>
#include <iostream>
+#include <math.h>
#include <numeric>
+#include <stdlib.h>
+#include <unistd.h>
+#include "FitEvent.h"
+#include "TGraph.h"
+#include "TH2Poly.h"
#include <TChain.h>
#include <TFile.h>
#include <TH1D.h>
#include <TH2D.h>
#include <THStack.h>
#include <TKey.h>
#include <TLegend.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TObjArray.h>
#include <TROOT.h>
#include <TRandom3.h>
#include <TTree.h>
-#include "FitEvent.h"
-#include "TGraph.h"
-#include "TH2Poly.h"
-#include "FitEvent.h"
#include "FitLogger.h"
/*!
* \addtogroup Utils
* @{
*/
/// Functions needed by individual samples for calculating kinematic quantities.
namespace FitUtils {
/// Return a vector of all values saved in map
double *GetArrayFromMap(std::vector<std::string> invals,
std::map<std::string, double> inmap);
/// Returns kinetic energy of particle
double T(TLorentzVector part);
/// Returns momentum of particle
double p(TLorentzVector part);
-double p(FitParticle* part);
+double p(FitParticle *part);
/// Returns angle between particles (_NOT_ cosine!)
double th(TLorentzVector part, TLorentzVector part2);
-double th(FitParticle* part1, FitParticle* part2);
+double th(FitParticle *part1, FitParticle *part2);
/// Hadronic mass reconstruction
double Wrec(TLorentzVector pnu, TLorentzVector pmu);
/// Hadronic mass true from initial state particles and muon; useful if the full
/// FSI vectors aren't not saved and we for some reasons need W_true
double Wtrue(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector pnuc);
double SumKE_PartVect(std::vector<FitParticle *> const fps);
double SumTE_PartVect(std::vector<FitParticle *> const fps);
/// Return E Hadronic for all FS Particles in Hadronic System
double GetErecoil_TRUE(FitEvent *event);
/// Return E Hadronic for all Charged FS Particles in Hadronic System
double GetErecoil_CHARGED(FitEvent *event);
double Eavailable(FitEvent *event);
/*
CCQE MiniBooNE/MINERvA
*/
/// Function to calculate the reconstructed Q^{2}_{QE}
-double Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true);
+double Q2QErec(TLorentzVector pmu, double costh, double binding,
+ bool neutrino = true);
/// Function returns the reconstructed E_{nu} values
-double EnuQErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true);
+double EnuQErec(TLorentzVector pmu, double costh, double binding,
+ bool neutrino = true);
/// Function returns the reconstructed E_{nu} values
-double EnuQErec(TLorentzVector pmu, TLorentzVector pnu, double binding, bool neutrino = true);
+double EnuQErec(TLorentzVector pmu, TLorentzVector pnu, double binding,
+ bool neutrino = true);
//! Function to calculate the reconstructed Q^{2}_{QE}
double Q2QErec(double pl, double costh, double binding, bool neutrino = true);
//! Function to calculate the reconstructed Q^{2}_{QE}
-double Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino = true);
+double Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding,
+ bool neutrino = true);
//! Function returns the reconstructed E_{nu} values
-double EnuQErec(double pl, double costh, double binding,
- bool neutrino = true);
+double EnuQErec(double pl, double costh, double binding, bool neutrino = true);
/*
CCQE1p MINERvA
*/
/// Reconstruct Q2QE given just the maximum energy proton.
double ProtonQ2QErec(double pE, double binding);
/*
E Recoil MINERvA
*/
double GetErecoil_MINERvA_LowRecoil(FitEvent *event);
/*
CC1pi0 MiniBooNE
*/
/// Reconstruct Enu from CCpi0 vectors and binding energy
double EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0));
/// Reconstruct Q2 from CCpi0 vectors and binding energy
double Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0));
/*
CC1pi+ MiniBooNE
*/
/// returns reconstructed Enu a la MiniBooNE CCpi+
/// returns reconstructed Enu a la MiniBooNE CCpi+
// Also for when not having pion info (so when we have a Michel tag in T2K)
double EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip,
bool pionInfo = true);
/// returns reconstructed Enu assumming resonance interaction where intermediate
/// resonance was a Delta
double EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu);
/// returns reconstructed in a variety of flavours
double Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip,
int enuType = 0, bool pionInfo = true);
/*
T2K CC1pi+ on CH
*/
double thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
double q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
double WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppip);
double EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
/*
nucleon single pion
*/
double MpPi(TLorentzVector pp, TLorentzVector ppi);
+const double kDefaultSTV_ProtonMinCut = 0;
+const double kDefaultSTV_ProtonMaxCut = 100E3;
+const double kDefaultSTV_PProtonCosThetaCut = -1;
+
/// Gets delta p T as defined in Phys.Rev. C94 (2016) no.1, 015503
-double Get_STV_dpt(FitEvent *event, double ProtonMinCut=0, double ProtonMaxCut=100E3, double ProtonCosThetaCut=-1, int ISPDG=14, bool Is0pi=true);
+double
+Get_STV_dpt_protonps(FitEvent *event,
+ double ProtonMinCut = kDefaultSTV_ProtonMinCut,
+ double ProtonMaxCut = kDefaultSTV_ProtonMaxCut,
+ double ProtonCosThetaCut = kDefaultSTV_PProtonCosThetaCut,
+ int ISPDG = 14, bool Is0pi = true);
/// Gets delta phi T as defined in Phys.Rev. C94 (2016) no.1, 015503
-double Get_STV_dphit(FitEvent *event, double ProtonMinCut=0, double ProtonMaxCut=100E3, double ProtonCosThetaCut=-1, int ISPDG=14, bool Is0pi=true);
+double Get_STV_dphit_protonps(
+ FitEvent *event, double ProtonMinCut = kDefaultSTV_ProtonMinCut,
+ double ProtonMaxCut = kDefaultSTV_ProtonMaxCut,
+ double ProtonCosThetaCut = kDefaultSTV_PProtonCosThetaCut, int ISPDG = 14,
+ bool Is0pi = true);
/// Gets delta alpha T as defined in Phys.Rev. C94 (2016) no.1, 015503
-double Get_STV_dalphat(FitEvent *event, double ProtonMinCut=0, double ProtonMaxCut=100E3, double ProtonCosThetaCut=-1, int ISPDG=14, bool Is0pi=true);
+double Get_STV_dalphat_protonps(
+ FitEvent *event, double ProtonMinCut = kDefaultSTV_ProtonMinCut,
+ double ProtonMaxCut = kDefaultSTV_ProtonMaxCut,
+ double ProtonCosThetaCut = kDefaultSTV_PProtonCosThetaCut, int ISPDG = 14,
+ bool Is0pi = true);
+
+inline double Get_STV_dpt(FitEvent *event, int ISPDG, bool Is0pi) {
+ return Get_STV_dpt_protonps(event, kDefaultSTV_ProtonMinCut,
+ kDefaultSTV_ProtonMaxCut,
+ kDefaultSTV_PProtonCosThetaCut, ISPDG, Is0pi);
+}
+inline double Get_STV_dphit(FitEvent *event, int ISPDG, bool Is0pi) {
+ return Get_STV_dphit_protonps(event, kDefaultSTV_ProtonMinCut,
+ kDefaultSTV_ProtonMaxCut,
+ kDefaultSTV_PProtonCosThetaCut, ISPDG, Is0pi);
+}
+inline double Get_STV_dalphat(FitEvent *event, int ISPDG, bool Is0pi) {
+ return Get_STV_dalphat_protonps(event, kDefaultSTV_ProtonMinCut,
+ kDefaultSTV_ProtonMaxCut,
+ kDefaultSTV_PProtonCosThetaCut, ISPDG, Is0pi);
+}
// As defined in PhysRevC.95.065501
-// Using prescription from arXiv 1805.05486
-double Get_pn_reco_C(FitEvent *event, double ProtonMinCut=0, double ProtonMaxCut=100E3, double ProtonCosThetaCut=-1, int ISPDG=14, bool Is0pi=true);
-double Get_pn_reco_Ar(FitEvent *event, double ProtonMinCut=0, double ProtonMaxCut=100E3, double ProtonCosThetaCut=-1, int ISPDG=14, bool Is0pi=true);
+// Using prescription from arXiv 1805.05486
+double Get_pn_reco_C_protonps(
+ FitEvent *event, double ProtonMinCut = kDefaultSTV_ProtonMinCut,
+ double ProtonMaxCut = kDefaultSTV_ProtonMaxCut,
+ double ProtonCosThetaCut = kDefaultSTV_PProtonCosThetaCut, int ISPDG = 14,
+ bool Is0pi = true);
+
+double Get_pn_reco_Ar_protonps(
+ FitEvent *event, double ProtonMinCut = kDefaultSTV_ProtonMinCut,
+ double ProtonMaxCut = kDefaultSTV_ProtonMaxCut,
+ double ProtonCosThetaCut = kDefaultSTV_PProtonCosThetaCut, int ISPDG = 14,
+ bool Is0pi = true);
+
+inline double Get_pn_reco_C(FitEvent *event, int ISPDG, bool Is0pi) {
+ return Get_pn_reco_C_protonps(event, kDefaultSTV_ProtonMinCut,
+ kDefaultSTV_ProtonMaxCut,
+ kDefaultSTV_PProtonCosThetaCut, ISPDG, Is0pi);
+}
-// Helper for STV; get a proton from an event within some cut range (different for T2K and MINERvA)
-TLorentzVector GetProtonInRange(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut);
+inline double Get_pn_reco_Ar(FitEvent *event, int ISPDG, bool Is0pi) {
+ return Get_pn_reco_Ar_protonps(event, kDefaultSTV_ProtonMinCut,
+ kDefaultSTV_ProtonMaxCut,
+ kDefaultSTV_PProtonCosThetaCut, ISPDG, Is0pi);
+}
-//For T2K inferred kinematics analyis - variables defined as on page 7 of T2K TN287v11 (and now arXiv 1802.05078)
-double ppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
-TVector3 tppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
-double cthpInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
+// Helper for STV; get a proton from an event within some cut range (different
+// for T2K and MINERvA)
+TLorentzVector GetProtonInRange(FitEvent *event, double ProtonMinCut,
+ double ProtonMaxCut, double ProtonCosThetaCut);
-double CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot);
-double PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot);
-}
+// For T2K inferred kinematics analyis - variables defined as on page 7 of T2K
+// TN287v11 (and now arXiv 1802.05078)
+double ppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
+TVector3 tppInfK(TLorentzVector pmu, double costh, double binding,
+ bool neutrino);
+double cthpInfK(TLorentzVector pmu, double costh, double binding,
+ bool neutrino);
+
+double CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi,
+ TLorentzVector Pprot);
+double PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi,
+ TLorentzVector Pprot);
+} // namespace FitUtils
/*! @} */
#endif