diff --git a/src/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
index 93adf7c..0e889ac 100644
--- a/src/MCStudies/GenericFlux_Vectors.cxx
+++ b/src/MCStudies/GenericFlux_Vectors.cxx
@@ -1,483 +1,483 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "GenericFlux_Vectors.h"
#ifndef __NO_MINERvA__
#include "MINERvA_SignalDef.h"
#endif
GenericFlux_Vectors::GenericFlux_Vectors(std::string name,
std::string inputfile, FitWeight *rw,
std::string type,
std::string fakeDataFile) {
// Measurement Details
fName = name;
eventVariables = NULL;
// Define our energy range for flux calcs
EnuMin = 0.;
EnuMax = 1E10; // Arbritrarily high energy limit
if (Config::HasPar("EnuMin")) {
EnuMin = Config::GetParD("EnuMin");
}
if (Config::HasPar("EnuMax")) {
EnuMax = Config::GetParD("EnuMax");
}
SavePreFSI = Config::Get().GetParB("nuisflat_SavePreFSI");
NUIS_LOG(SAM, "Running GenericFlux_Vectors saving pre-FSI particles? "
<< SavePreFSI);
// Set default fitter flags
fIsDiag = true;
fIsShape = false;
fIsRawEvents = false;
// This function will sort out the input files automatically and parse all the
// inputs,flags,etc.
// There may be complex cases where you have to do this by hand, but usually
// this will do.
Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
eventVariables = NULL;
// Setup fDataHist as a placeholder
this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
this->SetupDefaultHist();
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
covar = StatUtils::GetInvert(fFullCovar);
// 1. The generator is organised in SetupMeasurement so it gives the
// cross-section in "per nucleon" units.
// So some extra scaling for a specific measurement may be required. For
// Example to get a "per neutron" measurement on carbon
// which we do here, we have to multiple by the number of nucleons 12 and
// divide by the number of neutrons 6.
// N.B. MeasurementBase::PredictedEventRate includes the 1E-38 factor that is
// often included here in other classes that directly integrate the event
// histogram. This method is used here as it now respects EnuMin and EnuMax
// correctly.
this->fScaleFactor =
(this->PredictedEventRate("width", 0, EnuMax) / double(fNEvents)) /
this->TotalIntegratedFlux("width");
NUIS_LOG(SAM, " Generic Flux Scaling Factor = "
<< fScaleFactor << " [= "
<< (GetEventHistogram()->Integral("width") * 1E-38) << "/("
<< (fNEvents + 0.) << "*" << TotalIntegratedFlux("width")
<< ")]");
if (fScaleFactor <= 0.0) {
NUIS_ABORT("SCALE FACTOR TOO LOW");
}
// Setup our TTrees
this->AddEventVariablesToTree();
this->AddSignalFlagsToTree();
}
void GenericFlux_Vectors::AddEventVariablesToTree() {
// Setup the TTree to save everything
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
NUIS_LOG(SAM, "Adding Event Variables");
eventVariables->Branch("Mode", &Mode, "Mode/I");
eventVariables->Branch("cc", &cc, "cc/B");
eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I");
eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F");
eventVariables->Branch("tgt", &tgt, "tgt/I");
eventVariables->Branch("tgta", &tgta, "tgta/I");
eventVariables->Branch("tgtz", &tgtz, "tgtz/I");
eventVariables->Branch("PDGLep", &PDGLep, "PDGLep/I");
eventVariables->Branch("ELep", &ELep, "ELep/F");
eventVariables->Branch("CosLep", &CosLep, "CosLep/F");
// Basic interaction kinematics
eventVariables->Branch("Q2", &Q2, "Q2/F");
eventVariables->Branch("q0", &q0, "q0/F");
eventVariables->Branch("q3", &q3, "q3/F");
eventVariables->Branch("Enu_QE", &Enu_QE, "Enu_QE/F");
eventVariables->Branch("Q2_QE", &Q2_QE, "Q2_QE/F");
eventVariables->Branch("W_nuc_rest", &W_nuc_rest, "W_nuc_rest/F");
eventVariables->Branch("W", &W, "W/F");
eventVariables->Branch("W_genie", &W_genie, "W_genie/F");
eventVariables->Branch("x", &x, "x/F");
eventVariables->Branch("y", &y, "y/F");
eventVariables->Branch("Eav", &Eav, "Eav/F");
eventVariables->Branch("EavAlt", &EavAlt, "EavAlt/F");
eventVariables->Branch("CosThetaAdler", &CosThetaAdler, "CosThetaAdler/F");
eventVariables->Branch("PhiAdler", &PhiAdler, "PhiAdler/F");
eventVariables->Branch("dalphat", &dalphat, "dalphat/F");
eventVariables->Branch("dpt", &dpt, "dpt/F");
eventVariables->Branch("dphit", &dphit, "dphit/F");
eventVariables->Branch("pnreco_C", &pnreco_C, "pnreco_C/F");
// Save outgoing particle vectors
eventVariables->Branch("nfsp", &nfsp, "nfsp/I");
eventVariables->Branch("px", px, "px[nfsp]/F");
eventVariables->Branch("py", py, "py[nfsp]/F");
eventVariables->Branch("pz", pz, "pz[nfsp]/F");
eventVariables->Branch("E", E, "E[nfsp]/F");
eventVariables->Branch("pdg", pdg, "pdg[nfsp]/I");
eventVariables->Branch("pdg_rank", pdg_rank, "pdg_rank[nfsp]/I");
// Save init particle vectors
eventVariables->Branch("ninitp", &ninitp, "ninitp/I");
eventVariables->Branch("px_init", px_init, "px_init[ninitp]/F");
eventVariables->Branch("py_init", py_init, "py_init[ninitp]/F");
eventVariables->Branch("pz_init", pz_init, "pz_init[ninitp]/F");
eventVariables->Branch("E_init", E_init, "E_init[ninitp]/F");
eventVariables->Branch("pdg_init", pdg_init, "pdg_init[ninitp]/I");
// Save pre-FSI vectors
eventVariables->Branch("nvertp", &nvertp, "nvertp/I");
eventVariables->Branch("px_vert", px_vert, "px_vert[nvertp]/F");
eventVariables->Branch("py_vert", py_vert, "py_vert[nvertp]/F");
eventVariables->Branch("pz_vert", pz_vert, "pz_vert[nvertp]/F");
eventVariables->Branch("E_vert", E_vert, "E_vert[nvertp]/F");
eventVariables->Branch("pdg_vert", pdg_vert, "pdg_vert[nvertp]/I");
// Event Scaling Information
eventVariables->Branch("Weight", &Weight, "Weight/F");
eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
// Should be a double because may be 1E-39 and less
eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/D");
// The customs
eventVariables->Branch("CustomWeight", &CustomWeight, "CustomWeight/F");
eventVariables->Branch("CustomWeightArray", CustomWeightArray,
"CustomWeightArray[6]/F");
return;
}
void GenericFlux_Vectors::FillEventVariables(FitEvent *event) {
ResetVariables();
// Fill Signal Variables
FillSignalFlags(event);
NUIS_LOG(DEB, "Filling signal");
// Now fill the information
Mode = event->Mode;
cc = event->IsCC();
// Get the incoming neutrino and outgoing lepton
FitParticle *nu = event->GetBeamPart();
FitParticle *lep = event->GetHMFSAnyLepton();
PDGnu = nu->fPID;
Enu_true = nu->fP.E() / 1E3;
tgt = event->fTargetPDG;
tgta = event->fTargetA;
tgtz = event->fTargetZ;
TLorentzVector ISP4 = nu->fP;
if (lep != NULL) {
PDGLep = lep->fPID;
ELep = lep->fP.E() / 1E3;
CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect()));
// Basic interaction kinematics
Q2 = -1 * (nu->fP - lep->fP).Mag2() / 1E6;
q0 = (nu->fP - lep->fP).E() / 1E3;
q3 = (nu->fP - lep->fP).Vect().Mag() / 1E3;
// These assume C12 binding from MINERvA... not ideal
Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true);
Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true);
Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event) / 1.E3;
EavAlt = FitUtils::Eavailable(event) / 1.E3;
// Check if this is a 1pi+ or 1pi0 event
if ((SignalDef::isCC1pi(event, PDGnu, 211) ||
SignalDef::isCC1pi(event, PDGnu, -211) ||
SignalDef::isCC1pi(event, PDGnu, 111)) &&
event->NumFSNucleons() == 1) {
TLorentzVector Pnu = nu->fP;
TLorentzVector Pmu = lep->fP;
TLorentzVector Ppi = event->GetHMFSPions()->fP;
TLorentzVector Pprot = event->GetHMFSNucleons()->fP;
CosThetaAdler = FitUtils::CosThAdler(Pnu, Pmu, Ppi, Pprot);
PhiAdler = FitUtils::PhiAdler(Pnu, Pmu, Ppi, Pprot);
}
// Get W_true with assumption of initial state nucleon at rest
float m_n = (float)PhysConst::mass_proton;
// Q2 assuming nucleon at rest
W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
// True Q2
x = Q2 / (2 * m_n * q0);
y = 1 - ELep / Enu_true;
- dalphat = FitUtils::Get_STV_dalphat(event, PDGnu, true);
- dpt = FitUtils::Get_STV_dpt(event, PDGnu, true);
- dphit = FitUtils::Get_STV_dphit(event, PDGnu, true);
- pnreco_C = FitUtils::Get_pn_reco_C(event, PDGnu, true);
+ dalphat = FitUtils::Get_STV_dalphat_HMProton(event, PDGnu, true);
+ dpt = FitUtils::Get_STV_dpt_HMProton(event, PDGnu, true);
+ dphit = FitUtils::Get_STV_dphit_HMProton(event, PDGnu, true);
+ pnreco_C = FitUtils::Get_pn_reco_C_HMProton(event, PDGnu, true);
}
// Loop over the particles and store all the final state particles in a vector
for (UInt_t i = 0; i < event->Npart(); ++i) {
if (event->PartInfo(i)->fIsAlive &&
event->PartInfo(i)->Status() == kFinalState)
partList.push_back(event->PartInfo(i));
if (SavePreFSI && event->fPrimaryVertex[i])
vertList.push_back(event->PartInfo(i));
if (SavePreFSI && event->PartInfo(i)->IsInitialState())
initList.push_back(event->PartInfo(i));
if(event->PartInfo(i)->IsInitialState()){
ISP4 += event->PartInfo(i)->fP;
}
}
// Save outgoing particle vectors
nfsp = (int)partList.size();
std::map > > pdgMap;
for (int i = 0; i < nfsp; ++i) {
px[i] = partList[i]->fP.X() / 1E3;
py[i] = partList[i]->fP.Y() / 1E3;
pz[i] = partList[i]->fP.Z() / 1E3;
E[i] = partList[i]->fP.E() / 1E3;
pdg[i] = partList[i]->fPID;
pdgMap[pdg[i]].push_back(std::make_pair(partList[i]->fP.Vect().Mag(), i));
}
for (std::map > >::iterator iter =
pdgMap.begin();
iter != pdgMap.end(); ++iter) {
std::vector > thisVect = iter->second;
std::sort(thisVect.begin(), thisVect.end());
// Now save the order... a bit funky to avoid inverting
int nPart = (int)thisVect.size() - 1;
for (int i = nPart; i >= 0; --i) {
pdg_rank[thisVect[i].second] = nPart - i;
}
}
// Save pre-FSI particles
nvertp = (int)vertList.size();
for (int i = 0; i < nvertp; ++i) {
px_vert[i] = vertList[i]->fP.X() / 1E3;
py_vert[i] = vertList[i]->fP.Y() / 1E3;
pz_vert[i] = vertList[i]->fP.Z() / 1E3;
E_vert[i] = vertList[i]->fP.E() / 1E3;
pdg_vert[i] = vertList[i]->fPID;
}
// Save init particles
ninitp = (int)initList.size();
for (int i = 0; i < ninitp; ++i) {
px_init[i] = initList[i]->fP.X() / 1E3;
py_init[i] = initList[i]->fP.Y() / 1E3;
pz_init[i] = initList[i]->fP.Z() / 1E3;
E_init[i] = initList[i]->fP.E() / 1E3;
pdg_init[i] = initList[i]->fPID;
}
#ifdef __GENIE_ENABLED__
if (event->fType == kGENIE) {
EventRecord *gevent = static_cast(event->genie_event->event);
const Interaction *interaction = gevent->Summary();
const Kinematics &kine = interaction->Kine();
W_genie = kine.W();
}
#endif
if (lep != NULL) {
W = (ISP4 - lep->fP).M();
} else {
W = 0;
}
// Fill event weights
Weight = event->RWWeight * event->InputWeight;
RWWeight = event->RWWeight;
InputWeight = event->InputWeight;
// And the Customs
CustomWeight = event->CustomWeight;
for (int i = 0; i < 6; ++i) {
CustomWeightArray[i] = event->CustomWeightArray[i];
}
// Fill the eventVariables Tree
eventVariables->Fill();
return;
};
//********************************************************************
void GenericFlux_Vectors::ResetVariables() {
//********************************************************************
cc = false;
// Reset all Function used to extract any variables of interest to the event
Mode = PDGnu = tgt = tgta = tgtz = PDGLep = 0;
Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W =
x = y = Eav = EavAlt = CosThetaAdler = PhiAdler = -999.9;
W_genie = -999;
// Other fun variables
// MINERvA-like ones
dalphat = dpt = dphit = pnreco_C = -999.99;
nfsp = ninitp = nvertp = 0;
for (int i = 0; i < kMAX; ++i) {
px[i] = py[i] = pz[i] = E[i] = -999;
pdg[i] = pdg_rank[i] = 0;
px_init[i] = py_init[i] = pz_init[i] = E_init[i] = -999;
pdg_init[i] = 0;
px_vert[i] = py_vert[i] = pz_vert[i] = E_vert[i] = -999;
pdg_vert[i] = 0;
}
Weight = InputWeight = RWWeight = 0.0;
CustomWeight = 0.0;
for (int i = 0; i < 6; ++i)
CustomWeightArray[i] = 0.0;
partList.clear();
initList.clear();
vertList.clear();
flagCCINC = flagNCINC = flagCCQE = flagCC0pi = flagCCQELike = flagNCEL =
flagNC0pi = flagCCcoh = flagNCcoh = flagCC1pip = flagNC1pip = flagCC1pim =
flagNC1pim = flagCC1pi0 = flagNC1pi0 = false;
#ifndef __NO_MINERvA__
flagCC0piMINERvA = false;
#endif
}
//********************************************************************
void GenericFlux_Vectors::FillSignalFlags(FitEvent *event) {
//********************************************************************
// Some example flags are given from SignalDef.
// See src/Utils/SignalDef.cxx for more.
int nuPDG = event->PartInfo(0)->fPID;
// Generic signal flags
flagCCINC = SignalDef::isCCINC(event, nuPDG);
flagNCINC = SignalDef::isNCINC(event, nuPDG);
flagCCQE = SignalDef::isCCQE(event, nuPDG);
flagCCQELike = SignalDef::isCCQELike(event, nuPDG);
flagCC0pi = SignalDef::isCC0pi(event, nuPDG);
flagNCEL = SignalDef::isNCEL(event, nuPDG);
flagNC0pi = SignalDef::isNC0pi(event, nuPDG);
flagCCcoh = SignalDef::isCCCOH(event, nuPDG, 211);
flagNCcoh = SignalDef::isNCCOH(event, nuPDG, 111);
flagCC1pip = SignalDef::isCC1pi(event, nuPDG, 211);
flagNC1pip = SignalDef::isNC1pi(event, nuPDG, 211);
flagCC1pim = SignalDef::isCC1pi(event, nuPDG, -211);
flagNC1pim = SignalDef::isNC1pi(event, nuPDG, -211);
flagCC1pi0 = SignalDef::isCC1pi(event, nuPDG, 111);
flagNC1pi0 = SignalDef::isNC1pi(event, nuPDG, 111);
#ifndef __NO_MINERvA__
flagCC0piMINERvA = SignalDef::isCC0pi_MINERvAPTPZ(event, 14);
#endif
}
void GenericFlux_Vectors::AddSignalFlagsToTree() {
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
NUIS_LOG(SAM, "Adding signal flags");
// Signal Definitions from SignalDef.cxx
eventVariables->Branch("flagCCINC", &flagCCINC, "flagCCINC/O");
eventVariables->Branch("flagNCINC", &flagNCINC, "flagNCINC/O");
eventVariables->Branch("flagCCQE", &flagCCQE, "flagCCQE/O");
eventVariables->Branch("flagCC0pi", &flagCC0pi, "flagCC0pi/O");
eventVariables->Branch("flagCCQELike", &flagCCQELike, "flagCCQELike/O");
eventVariables->Branch("flagNCEL", &flagNCEL, "flagNCEL/O");
eventVariables->Branch("flagNC0pi", &flagNC0pi, "flagNC0pi/O");
eventVariables->Branch("flagCCcoh", &flagCCcoh, "flagCCcoh/O");
eventVariables->Branch("flagNCcoh", &flagNCcoh, "flagNCcoh/O");
eventVariables->Branch("flagCC1pip", &flagCC1pip, "flagCC1pip/O");
eventVariables->Branch("flagNC1pip", &flagNC1pip, "flagNC1pip/O");
eventVariables->Branch("flagCC1pim", &flagCC1pim, "flagCC1pim/O");
eventVariables->Branch("flagNC1pim", &flagNC1pim, "flagNC1pim/O");
eventVariables->Branch("flagCC1pi0", &flagCC1pi0, "flagCC1pi0/O");
eventVariables->Branch("flagNC1pi0", &flagNC1pi0, "flagNC1pi0/O");
#ifndef __NO_MINERvA__
eventVariables->Branch("flagCC0piMINERvA", &flagCC0piMINERvA,
"flagCC0piMINERvA/O");
#endif
};
void GenericFlux_Vectors::Write(std::string drawOpt) {
// First save the TTree
eventVariables->Write();
// Save Flux and Event Histograms too
GetInput()->GetFluxHistogram()->Write();
GetInput()->GetEventHistogram()->Write();
return;
}
// Override functions which aren't really necessary
bool GenericFlux_Vectors::isSignal(FitEvent *event) {
(void)event;
return true;
};
void GenericFlux_Vectors::ScaleEvents() { return; }
void GenericFlux_Vectors::ApplyNormScale(float norm) {
this->fCurrentNorm = norm;
return;
}
void GenericFlux_Vectors::FillHistograms() { return; }
void GenericFlux_Vectors::ResetAll() {
eventVariables->Reset();
return;
}
float GenericFlux_Vectors::GetChi2() { return 0.0; }
diff --git a/src/MINERvA/MINERvAUtils.cxx b/src/MINERvA/MINERvAUtils.cxx
index c9d20e2..63cf721 100644
--- a/src/MINERvA/MINERvAUtils.cxx
+++ b/src/MINERvA/MINERvAUtils.cxx
@@ -1,312 +1,571 @@
// 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 .
-*******************************************************************************/
+ * This file is part of NUISANCE.
+ *
+ * NUISANCE is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * NUISANCE is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with NUISANCE. If not, see .
+ *******************************************************************************/
#include "MINERvAUtils.h"
#include "FitUtils.h"
-namespace MINERvAPar{
- double MINERvADensity = FitPar::Config().GetParD("MINERvADensity");
- double MINERvARecoDist = FitPar::Config().GetParD("MINERvARecoDist");
- double PenetratingMuonE = FitPar::Config().GetParD("PenetratingMuonEnergy");
- double NumRangeSteps = FitPar::Config().GetParI("NumRangeSteps");
-}
+namespace MINERvAPar {
+double MINERvADensity = FitPar::Config().GetParD("MINERvADensity");
+double MINERvARecoDist = FitPar::Config().GetParD("MINERvARecoDist");
+double PenetratingMuonE = FitPar::Config().GetParD("PenetratingMuonEnergy");
+double NumRangeSteps = FitPar::Config().GetParI("NumRangeSteps");
+} // namespace MINERvAPar
-double MINERvAUtils::StoppedEfficiency(TH2D *effHist, FitParticle *nu, FitParticle *muon){
+double MINERvAUtils::StoppedEfficiency(TH2D *effHist, FitParticle *nu,
+ FitParticle *muon) {
double eff = 0.;
- if (!effHist) return eff;
- eff = effHist->GetBinContent(effHist->FindBin(FitUtils::p(muon), FitUtils::th(nu, muon)/TMath::Pi()*180.));
+ if (!effHist)
+ return eff;
+ eff = effHist->GetBinContent(effHist->FindBin(
+ FitUtils::p(muon), FitUtils::th(nu, muon) / TMath::Pi() * 180.));
return eff;
}
-double MINERvAUtils::PenetratedEfficiency(FitParticle *nu, FitParticle *muon){
+double MINERvAUtils::PenetratedEfficiency(FitParticle *nu, FitParticle *muon) {
double eff = 0.;
- if (FitUtils::th(nu, muon)/TMath::Pi()*180. > 50) eff = 0.;
- if (FitUtils::p(muon) < 1.4) eff = 0.;
+ if (FitUtils::th(nu, muon) / TMath::Pi() * 180. > 50)
+ eff = 0.;
+ if (FitUtils::p(muon) < 1.4)
+ eff = 0.;
return eff;
}
-double MINERvAUtils::BetheBlochCH(double E, double mass){
+double MINERvAUtils::BetheBlochCH(double E, double mass) {
- double beta2 = 1 - mass*mass/E/E;
- double gamma = 1./sqrt(1-beta2);
- double mass_ratio = PhysConst::mass_electron*1000./mass;
+ double beta2 = 1 - mass * mass / E / E;
+ double gamma = 1. / sqrt(1 - beta2);
+ double mass_ratio = PhysConst::mass_electron * 1000. / mass;
// Argh, have to remember to convert to MeV or you'll hate yourself!
- double I2 = 68.7e-6*68.7e-6; // mean excitation potential I
-
-
- double w_max = 2*PhysConst::mass_electron*1000.*beta2*gamma*gamma;
- w_max /= 1 + 2*gamma*mass_ratio + mass_ratio*mass_ratio;
+ double I2 = 68.7e-6 * 68.7e-6; // mean excitation potential I
+ double w_max = 2 * PhysConst::mass_electron * 1000. * beta2 * gamma * gamma;
+ w_max /= 1 + 2 * gamma * mass_ratio + mass_ratio * mass_ratio;
- // Values taken from the PDG for K = 0.307075 MeV mol-1 cm2, mean ionization energy I = 68.7 eV (Polystyrene)
- // = 0.53768 (pdg.lbl.gov/AtomicNuclearProperties)
- double log_term = log(2*PhysConst::mass_electron*1000.*beta2*gamma*gamma*w_max/I2);
- double dedx = 0.307075*0.53768/beta2*(0.5*log_term - beta2);
+ // Values taken from the PDG for K = 0.307075 MeV mol-1 cm2, mean ionization
+ // energy I = 68.7 eV (Polystyrene) = 0.53768
+ // (pdg.lbl.gov/AtomicNuclearProperties)
+ double log_term = log(2 * PhysConst::mass_electron * 1000. * beta2 * gamma *
+ gamma * w_max / I2);
+ double dedx = 0.307075 * 0.53768 / beta2 * (0.5 * log_term - beta2);
return dedx;
}
-
-// This function returns an estimate of the range of the particle in scintillator.
-// It uses crude integration and Bethe-Bloch to approximate the range.
-double MINERvAUtils::RangeInScintillator(FitParticle* particle, int nsteps){
+// This function returns an estimate of the range of the particle in
+// scintillator. It uses crude integration and Bethe-Bloch to approximate the
+// range.
+double MINERvAUtils::RangeInScintillator(FitParticle *particle, int nsteps) {
// The particle energy
- double E = particle->fP.E();
- double M = particle->fP.M();
+ double E = particle->fP.E();
+ double M = particle->fP.M();
double Ek = E - M;
- double step_size = Ek/float(nsteps+1);
+ double step_size = Ek / float(nsteps + 1);
double range = 0;
// Add an offset to make the integral a touch more accurate
- Ek -= step_size/2.;
- for (int i = 0; i < nsteps; ++i){
-
- double dEdx = MINERvAUtils::BetheBlochCH(Ek+M, M);
-
+ Ek -= step_size / 2.;
+ for (int i = 0; i < nsteps; ++i) {
+
+ double dEdx = MINERvAUtils::BetheBlochCH(Ek + M, M);
+
Ek -= step_size;
// dEdx is -ve
- range -= step_size/dEdx;
+ range -= step_size / dEdx;
}
// Account for density of polystyrene
range /= MINERvAPar::MINERvADensity;
// Range estimate is in cm
return fabs(range);
}
-double MINERvAUtils::GetEDepositOutsideRangeInScintillator(FitParticle* particle, double rangelimit){
+double
+MINERvAUtils::GetEDepositOutsideRangeInScintillator(FitParticle *particle,
+ double rangelimit) {
// The particle energy
- double E = particle->fP.E();
- double M = particle->fP.M();
+ double E = particle->fP.E();
+ double M = particle->fP.M();
double Ek = E - M;
int nsteps = MINERvAPar::NumRangeSteps;
- double step_size = Ek/float(nsteps+1);
+ double step_size = Ek / float(nsteps + 1);
double range = 0;
double Ekinside = 0.0;
double Ekstart = Ek;
// Add an offset to make the integral a touch more accurate
- Ek -= step_size/2.;
- for (int i = 0; i < nsteps; ++i){
-
- double dEdx = MINERvAUtils::BetheBlochCH(Ek+M, M);
-
+ Ek -= step_size / 2.;
+ for (int i = 0; i < nsteps; ++i) {
+
+ double dEdx = MINERvAUtils::BetheBlochCH(Ek + M, M);
+
Ek -= step_size;
// dEdx is -ve
- range -= step_size/dEdx;
- if (fabs(range) / MINERvAPar::MINERvADensity < rangelimit){
+ range -= step_size / dEdx;
+ if (fabs(range) / MINERvAPar::MINERvADensity < rangelimit) {
Ekinside = Ek;
}
}
// Range estimate is in cm
return Ekstart - Ekinside;
-}
+}
-double MINERvAUtils::GetEDepositInsideRangeInScintillator(FitParticle* particle, double rangelimit){
+double MINERvAUtils::GetEDepositInsideRangeInScintillator(FitParticle *particle,
+ double rangelimit) {
// The particle energy
- double E = particle->fP.E();
- double M = particle->fP.M();
+ double E = particle->fP.E();
+ double M = particle->fP.M();
double Ek = E - M;
return Ek - GetEDepositOutsideRangeInScintillator(particle, rangelimit);
}
-
-
// Function to calculate the distance the particle travels in scintillator
-bool MINERvAUtils::PassesDistanceCut(FitParticle* beam, FitParticle* particle){
+bool MINERvAUtils::PassesDistanceCut(FitParticle *beam, FitParticle *particle) {
- double dist = MINERvAUtils::RangeInScintillator(particle, MINERvAPar::NumRangeSteps);
- double zdist = dist*cos(FitUtils::th(beam, particle));
+ double dist =
+ MINERvAUtils::RangeInScintillator(particle, MINERvAPar::NumRangeSteps);
+ double zdist = dist * cos(FitUtils::th(beam, particle));
- if (abs(zdist) < MINERvAPar::MINERvARecoDist) return false;
+ if (abs(zdist) < MINERvAPar::MINERvARecoDist)
+ return false;
return true;
}
-
// Function to return the MainTrk
-int MINERvAUtils::GetMainTrack(FitEvent *event, TH2D *effHist, FitParticle*& mainTrk, double& weight, bool penetrated){
+int MINERvAUtils::GetMainTrack(FitEvent *event, TH2D *effHist,
+ FitParticle *&mainTrk, double &weight,
+ bool penetrated) {
- FitParticle *nu = event->GetNeutrinoIn();
+ FitParticle *nu = event->GetNeutrinoIn();
double highestMom = 0;
int index = 0;
mainTrk = NULL;
// Loop over particles
- for (uint j = 2; j < event->Npart(); ++j){
+ for (uint j = 2; j < event->Npart(); ++j) {
// Final state only!
- if (!(event->PartInfo(j))->fIsAlive) continue;
- if (event->PartInfo(j)->fNEUTStatusCode != 0) continue;
+ if (!(event->PartInfo(j))->fIsAlive)
+ continue;
+ if (event->PartInfo(j)->fNEUTStatusCode != 0)
+ continue;
int PID = event->PartInfo(j)->fPID;
// Only consider pions, muons for now
- if (abs(PID) != 211 && abs(PID) != 13) continue;
+ if (abs(PID) != 211 && abs(PID) != 13)
+ continue;
// Ignore if higher momentum tracks available
- if (FitUtils::p(event->PartInfo(j)) < highestMom) continue;
+ if (FitUtils::p(event->PartInfo(j)) < highestMom)
+ continue;
// Okay, now this is highest momentum
highestMom = FitUtils::p(event->PartInfo(j));
- weight *= MINERvAUtils::StoppedEfficiency(effHist, nu, event->PartInfo(j));
- index = j;
+ weight *= MINERvAUtils::StoppedEfficiency(effHist, nu, event->PartInfo(j));
+ index = j;
mainTrk = event->PartInfo(j);
} // end loop over particle stack
-
+
return index;
}
-
-void MINERvAUtils::GetOtherTrackInfo(FitEvent *event, int mainIndex, int& nProtons, int& nPiMus, int& nVertex, FitParticle*& secondTrk){
+void MINERvAUtils::GetOtherTrackInfo(FitEvent *event, int mainIndex,
+ int &nProtons, int &nPiMus, int &nVertex,
+ FitParticle *&secondTrk) {
// Reset everything
- nPiMus = 0;
- nProtons = 0;
- nVertex = 0;
- secondTrk = NULL;
+ nPiMus = 0;
+ nProtons = 0;
+ nVertex = 0;
+ secondTrk = NULL;
- double highestMom = 0.;
+ double highestMom = 0.;
// Loop over particles
- for (uint j = 2; j < event->Npart(); ++j){
+ for (uint j = 2; j < event->Npart(); ++j) {
// Don't re-count the main track
- if (j == (uint)mainIndex) continue;
+ if (j == (uint)mainIndex)
+ continue;
// Final state only!
- if (!(event->PartInfo(j))->fIsAlive) continue;
- if (event->PartInfo(j)->fNEUTStatusCode != 0) continue;
+ if (!(event->PartInfo(j))->fIsAlive)
+ continue;
+ if (event->PartInfo(j)->fNEUTStatusCode != 0)
+ continue;
int PID = event->PartInfo(j)->fPID;
// Only consider pions, muons, protons
- if (abs(PID) != 211 && PID != 2212 && abs(PID) != 13) continue;
+ if (abs(PID) != 211 && PID != 2212 && abs(PID) != 13)
+ continue;
// Must be reconstructed as a track in SciBooNE
- if (MINERvAUtils::PassesDistanceCut(event->PartInfo(0), event->PartInfo(j))){
+ if (MINERvAUtils::PassesDistanceCut(event->PartInfo(0),
+ event->PartInfo(j))) {
// Keep track of the second highest momentum track
- if (FitUtils::p(event->PartInfo(j)) > highestMom){
- highestMom = FitUtils::p(event->PartInfo(j));
- secondTrk = event->PartInfo(j);
+ if (FitUtils::p(event->PartInfo(j)) > highestMom) {
+ highestMom = FitUtils::p(event->PartInfo(j));
+ secondTrk = event->PartInfo(j);
}
- if (PID == 2212) nProtons += 1;
- else nPiMus += 1;
- } else nVertex += 1;
+ if (PID == 2212)
+ nProtons += 1;
+ else
+ nPiMus += 1;
+ } else
+ nVertex += 1;
} // end loop over particle stack
return;
}
-
// NOTE: need to adapt this to allow for penetrating events...
// Simpler, but gives the same results as in Hirade-san's thesis
-double MINERvAUtils::CalcThetaPr(FitEvent *event, FitParticle *main, FitParticle *second, bool penetrated){
-
- FitParticle *nu = event->GetNeutrinoIn();
+double MINERvAUtils::CalcThetaPr(FitEvent *event, FitParticle *main,
+ FitParticle *second, bool penetrated) {
+
+ FitParticle *nu = event->GetNeutrinoIn();
- if (!main || !nu || !second) return -999;
+ if (!main || !nu || !second)
+ return -999;
// Construct the vector p_pr = (-p_mux, -p_muy, Enurec - pmucosthetamu)
- // where p_mux, p_muy are the projections of the candidate muon momentum onto the x and y dimension respectively
- double pmu = main->fP.Vect().Mag();
+ // where p_mux, p_muy are the projections of the candidate muon momentum onto
+ // the x and y dimension respectively
+ double pmu = main->fP.Vect().Mag();
double pmu_x = main->fP.Vect().X();
double pmu_y = main->fP.Vect().Y();
- if (penetrated){
+ if (penetrated) {
pmu = 1400.;
- double ratio = 1.4/main->fP.Vect().Mag();
- TVector3 mod_mu = main->fP.Vect()*ratio;
+ double ratio = 1.4 / main->fP.Vect().Mag();
+ TVector3 mod_mu = main->fP.Vect() * ratio;
pmu_x = mod_mu.X();
pmu_y = mod_mu.Y();
}
- double Enuqe = FitUtils::EnuQErec(pmu/1000.,cos(FitUtils::th(nu, main)), 27., true)*1000.;
- double p_pr_z = Enuqe - pmu*cos(FitUtils::th(nu, main));
+ double Enuqe =
+ FitUtils::EnuQErec(pmu / 1000., cos(FitUtils::th(nu, main)), 27., true) *
+ 1000.;
+ double p_pr_z = Enuqe - pmu * cos(FitUtils::th(nu, main));
- TVector3 p_pr = TVector3(-pmu_x, -pmu_y, p_pr_z);
- double thetapr = p_pr.Angle(second->fP.Vect())/TMath::Pi()*180.;
+ TVector3 p_pr = TVector3(-pmu_x, -pmu_y, p_pr_z);
+ double thetapr = p_pr.Angle(second->fP.Vect()) / TMath::Pi() * 180.;
return thetapr;
}
-double MINERvAUtils::CalcThetaPi(FitEvent *event, FitParticle *second){
+double MINERvAUtils::CalcThetaPi(FitEvent *event, FitParticle *second) {
- FitParticle *nu = event->GetNeutrinoIn();
+ FitParticle *nu = event->GetNeutrinoIn();
- if (!second || !nu) return -999;
+ if (!second || !nu)
+ return -999;
- double thetapi = FitUtils::th(nu, second)/TMath::Pi()*180.;
+ double thetapi = FitUtils::th(nu, second) / TMath::Pi() * 180.;
return thetapi;
}
/// Functions to deal with the SB mode stacks
-MINERvAUtils::ModeStack::ModeStack(std::string name, std::string title, TH1* hist) {
+MINERvAUtils::ModeStack::ModeStack(std::string name, std::string title,
+ TH1 *hist) {
fName = name;
fTitle = title;
- AddMode(0, "Other", "Other", kGreen+2, 2, 3244);
- AddMode(1, "CCDIS", "#nu_{#mu} CC DIS", kRed, 2, 3304);
- AddMode(2, "CCRES", "#nu_{#mu} CC Resonant", kGray+2, 2, 1001);
- AddMode(3, "CCQE", "#nu_{#mu} CC QE", kMagenta, 2, 1001);
- AddMode(4, "CC2P2H", "#nu_{#mu} CC 2p2h", kAzure+1, 2, 1001);
+ AddMode(0, "Other", "Other", kGreen + 2, 2, 3244);
+ AddMode(1, "CCDIS", "#nu_{#mu} CC DIS", kRed, 2, 3304);
+ AddMode(2, "CCRES", "#nu_{#mu} CC Resonant", kGray + 2, 2, 1001);
+ AddMode(3, "CCQE", "#nu_{#mu} CC QE", kMagenta, 2, 1001);
+ AddMode(4, "CC2P2H", "#nu_{#mu} CC 2p2h", kAzure + 1, 2, 1001);
StackBase::SetupStack(hist);
};
-int MINERvAUtils::ModeStack::ConvertModeToIndex(int mode){
- switch (abs(mode)){
- case 1: return 3;
- case 2: return 4;
+int MINERvAUtils::ModeStack::ConvertModeToIndex(int mode) {
+ switch (abs(mode)) {
+ case 1:
+ return 3;
+ case 2:
+ return 4;
case 11:
case 12:
- case 13: return 2;
- case 26: return 1;
- default: return 0;
+ case 13:
+ return 2;
+ case 26:
+ return 1;
+ default:
+ return 0;
}
};
-
-void MINERvAUtils::ModeStack::Fill(int mode, double x, double y, double z, double weight) {
- StackBase::FillStack(MINERvAUtils::ModeStack::ConvertModeToIndex(mode), x, y, z, weight);
+void MINERvAUtils::ModeStack::Fill(int mode, double x, double y, double z,
+ double weight) {
+ StackBase::FillStack(MINERvAUtils::ModeStack::ConvertModeToIndex(mode), x, y,
+ z, weight);
};
-void MINERvAUtils::ModeStack::Fill(FitEvent* evt, double x, double y, double z, double weight) {
- StackBase::FillStack(MINERvAUtils::ModeStack::ConvertModeToIndex(evt->Mode), x, y, z, weight);
+void MINERvAUtils::ModeStack::Fill(FitEvent *evt, double x, double y, double z,
+ double weight) {
+ StackBase::FillStack(MINERvAUtils::ModeStack::ConvertModeToIndex(evt->Mode),
+ x, y, z, weight);
};
-void MINERvAUtils::ModeStack::Fill(BaseFitEvt* evt, double x, double y, double z, double weight) {
- StackBase::FillStack(MINERvAUtils::ModeStack::ConvertModeToIndex(evt->Mode), x, y, z, weight);
+void MINERvAUtils::ModeStack::Fill(BaseFitEvt *evt, double x, double y,
+ double z, double weight) {
+ StackBase::FillStack(MINERvAUtils::ModeStack::ConvertModeToIndex(evt->Mode),
+ x, y, z, weight);
};
+
+double MINERvAUtils::Get_STV_dpt_MINERvAPS(FitEvent *event,
+ double ProtonMinCut_MeV,
+ double ProtonMaxCut_MeV,
+ 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_MeV and
+ // ProtonMaxCut_MeV MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = MINERvAUtils::GetProtonInRange(
+ event, ProtonMinCut_MeV, ProtonMaxCut_MeV, 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 FitUtils::GetDeltaPT(LeptonP, HadronP, NuP).Mag();
+}
+
+double MINERvAUtils::Get_STV_dphit_MINERvAPS(FitEvent *event,
+ double ProtonMinCut_MeV,
+ double ProtonMaxCut_MeV,
+ 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_MeV and
+ // ProtonMaxCut_MeV MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = MINERvAUtils::GetProtonInRange(
+ event, ProtonMinCut_MeV, ProtonMaxCut_MeV, 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 FitUtils::GetDeltaPhiT(LeptonP, HadronP, NuP);
+}
+
+double MINERvAUtils::Get_STV_dalphat_MINERvAPS(FitEvent *event,
+ double ProtonMinCut_MeV,
+ double ProtonMaxCut_MeV,
+ 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_MeV and
+ // ProtonMaxCut_MeV MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = MINERvAUtils::GetProtonInRange(
+ event, ProtonMinCut_MeV, ProtonMaxCut_MeV, 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 FitUtils::GetDeltaAlphaT(LeptonP, HadronP, NuP);
+}
+
+// As defined in PhysRevC.95.065501
+// Using prescription from arXiv 1805.05486
+// Returns in GeV
+double MINERvAUtils::Get_pn_reco_C_MINERvAPS(FitEvent *event,
+ double ProtonMinCut_MeV,
+ double ProtonMaxCut_MeV,
+ double ProtonCosThetaCut,
+ int ISPDG, bool Is0pi) {
+
+ 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_MeV and
+ // ProtonMaxCut_MeV MeV with cos(theta_p) < ProtonCosThetaCut
+ TLorentzVector Pprot = MINERvAUtils::GetProtonInRange(
+ event, ProtonMinCut_MeV, ProtonMaxCut_MeV, ProtonCosThetaCut);
+ TVector3 HadronP = Pprot.Vect();
+
+ 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 = FitUtils::GetDeltaPT(LeptonP, HadronP, NuP);
+ 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 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 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;
+
+ return pn_reco;
+}
+
+// Find the highest momentum proton in the event between ProtonMinCut_MeV and
+// ProtonMaxCut_MeV MeV with cos(theta_p) > ProtonCosThetaCut
+TLorentzVector MINERvAUtils::GetProtonInRange(FitEvent *event,
+ double ProtonMinCut_MeV,
+ double ProtonMaxCut_MeV,
+ double ProtonCosThetaCut) {
+ // Get the neutrino
+ TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
+ int HMFSProton = -1;
+ double HighestMomentum = 0.0;
+ // Get the stack of protons
+ std::vector Protons = event->GetAllFSProton();
+ for (size_t i = 0; i < Protons.size(); ++i) {
+ if (Protons[i]->p() > ProtonMinCut_MeV &&
+ Protons[i]->p() < ProtonMaxCut_MeV &&
+ cos(Protons[i]->P3().Angle(Pnu.Vect())) > ProtonCosThetaCut &&
+ Protons[i]->p() > HighestMomentum) {
+ HighestMomentum = Protons[i]->p();
+ HMFSProton = i;
+ }
+ }
+ if (HMFSProton == -1) {
+ return TLorentzVector(0, 0, 0, 0);
+ }
+ // Now get the proton
+ TLorentzVector Pprot = Protons[HMFSProton]->fP;
+ return Pprot;
+}
diff --git a/src/MINERvA/MINERvAUtils.h b/src/MINERvA/MINERvAUtils.h
index 890a291..5aa3b25 100644
--- a/src/MINERvA/MINERvAUtils.h
+++ b/src/MINERvA/MINERvAUtils.h
@@ -1,107 +1,139 @@
// 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 .
-*******************************************************************************/
+ * This file is part of NUISANCE.
+ *
+ * NUISANCE is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * NUISANCE is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with NUISANCE. If not, see .
+ *******************************************************************************/
#ifndef MINERvAUtils_H_SEEN
#define MINERvAUtils_H_SEEN
-#include
-#include
-#include
+#include
+#include
#include
#include
#include
-#include
-#include
#include
+#include
+#include
#include
+#include
#include
+#include
#include
-#include
-#include
-#include "FitParticle.h"
#include "FitEvent.h"
#include "FitLogger.h"
+#include "FitParticle.h"
#include "StandardStacks.h"
/*!
* \addtogroup Utils
* @{
*/
namespace MINERvAPar {
- extern double MINERvADensity;
- extern double MINERvARecoDist;
- extern double PenetratingMuonE;
- extern double NumRangeSteps;
-}
+extern double MINERvADensity;
+extern double MINERvARecoDist;
+extern double PenetratingMuonE;
+extern double NumRangeSteps;
+} // namespace MINERvAPar
namespace MINERvAUtils {
- double StoppedEfficiency(TH2D *effHist, FitParticle *nu, FitParticle *muon);
- double PenetratedEfficiency(FitParticle *nu, FitParticle *muon);
- double BetheBlochCH(double beta, double mass);
- double RangeInScintillator(FitParticle* particle, int nsteps=50);
- double GetEDepositOutsideRangeInScintillator(FitParticle* particle, double rangelimit);
- double GetEDepositInsideRangeInScintillator(FitParticle* particle, double rangelimit);
- bool PassesDistanceCut(FitParticle* beam, FitParticle* particle);
-
- int GetMainTrack(FitEvent *event, TH2D *effHist, FitParticle*& mainTrk, double& weight, bool penetrated=false);
- void GetOtherTrackInfo(FitEvent *event, int mainIndex, int& nProtons, int& nPiMus, int& nVertex, FitParticle*& secondTrk);
-
- double CalcThetaPr(FitEvent *event, FitParticle *main, FitParticle *second, bool penetrated=false);
- double CalcThetaPi(FitEvent *event, FitParticle *second);
-
- /// Break down the plots as in the MINERvA CCQE Papers
- class ModeStack : public StackBase {
- public:
-
- /// Main constructor listing true mode categories.
- ModeStack(std::string name, std::string title, TH1* hist);
-
- /// List to convert Modes to Index.
- /// Should be kept in sync with constructor.
- int ConvertModeToIndex(int mode);
- /// Fill from given mode integer
- void Fill(int mode, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
- /// Extracts Mode from FitEvent and fills
- void Fill(FitEvent* evt, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
- /// Extracts Mode from BaseFitEvt
- void Fill(BaseFitEvt* evt, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
-
- };
-
- /// Break down the plots in terms of 1-N pion contributions
- class PionStack : public StackBase {
- public:
-
- /// Main constructor listing true mode categories.
- PionStack(std::string name, std::string title, TH1* hist);
-
- /// List to convert Pions to Index.
- /// Should be kept in sync with constructor.
- int ConvertNPionsToIndex(int npions);
-
- /// Fill from given mode integer
- void Fill(int npions, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
- };
-
-
-}
+double StoppedEfficiency(TH2D *effHist, FitParticle *nu, FitParticle *muon);
+double PenetratedEfficiency(FitParticle *nu, FitParticle *muon);
+double BetheBlochCH(double beta, double mass);
+double RangeInScintillator(FitParticle *particle, int nsteps = 50);
+double GetEDepositOutsideRangeInScintillator(FitParticle *particle,
+ double rangelimit);
+double GetEDepositInsideRangeInScintillator(FitParticle *particle,
+ double rangelimit);
+bool PassesDistanceCut(FitParticle *beam, FitParticle *particle);
+
+int GetMainTrack(FitEvent *event, TH2D *effHist, FitParticle *&mainTrk,
+ double &weight, bool penetrated = false);
+void GetOtherTrackInfo(FitEvent *event, int mainIndex, int &nProtons,
+ int &nPiMus, int &nVertex, FitParticle *&secondTrk);
+
+double CalcThetaPr(FitEvent *event, FitParticle *main, FitParticle *second,
+ bool penetrated = false);
+double CalcThetaPi(FitEvent *event, FitParticle *second);
+
+/// Gets delta p T as defined in Phys.Rev. C94 (2016) no.1, 015503
+double Get_STV_dpt_MINERvAPS(FitEvent *event, double ProtonMinCut_MeV,
+ double ProtonMaxCut_MeV, double ProtonCosThetaCut,
+ int ISPDG = 14, bool Is0pi = true);
+/// Gets delta phi T as defined in Phys.Rev. C94 (2016) no.1, 015503
+double Get_STV_dphit_MINERvAPS(FitEvent *event, double ProtonMinCut_MeV,
+ double ProtonMaxCut_MeV,
+ double ProtonCosThetaCut, int ISPDG = 14,
+ bool Is0pi = true);
+/// Gets delta alpha T as defined in Phys.Rev. C94 (2016) no.1, 015503
+double Get_STV_dalphat_MINERvAPS(FitEvent *event, double ProtonMinCut_MeV,
+ double ProtonMaxCut_MeV,
+ double ProtonCosThetaCut, int ISPDG = 14,
+ bool Is0pi = true);
+
+// As defined in PhysRevC.95.065501
+// Using prescription from arXiv 1805.05486
+double Get_pn_reco_C_MINERvAPS(FitEvent *event, double ProtonMinCut_MeV,
+ double ProtonMaxCut_MeV,
+ double ProtonCosThetaCut, int ISPDG = 14,
+ bool Is0pi = true);
+
+// Helper for STV; get a proton from an event within some cut range (different
+// for T2K and MINERvA)
+TLorentzVector GetProtonInRange(FitEvent *event, double ProtonMinCut_MeV,
+ double ProtonMaxCut_MeV,
+ double ProtonCosThetaCut);
+
+/// Break down the plots as in the MINERvA CCQE Papers
+class ModeStack : public StackBase {
+public:
+ /// Main constructor listing true mode categories.
+ ModeStack(std::string name, std::string title, TH1 *hist);
+
+ /// List to convert Modes to Index.
+ /// Should be kept in sync with constructor.
+ int ConvertModeToIndex(int mode);
+ /// Fill from given mode integer
+ void Fill(int mode, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
+ /// Extracts Mode from FitEvent and fills
+ void Fill(FitEvent *evt, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
+ /// Extracts Mode from BaseFitEvt
+ void Fill(BaseFitEvt *evt, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
+};
+
+/// Break down the plots in terms of 1-N pion contributions
+class PionStack : public StackBase {
+public:
+ /// Main constructor listing true mode categories.
+ PionStack(std::string name, std::string title, TH1 *hist);
+
+ /// List to convert Pions to Index.
+ /// Should be kept in sync with constructor.
+ int ConvertNPionsToIndex(int npions);
+
+ /// Fill from given mode integer
+ void Fill(int npions, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
+};
+
+} // namespace MINERvAUtils
#endif
diff --git a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx
index b1b34cb..25117fe 100644
--- a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx
+++ b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx
@@ -1,353 +1,354 @@
// 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 .
*******************************************************************************/
// 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"
+#include "MINERvAUtils.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 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 CrossSection;
std::vector Error;
std::vector 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 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));
+ TLorentzVector Pprot = MINERvAUtils::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_protonps(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true);
+ fXVar = MINERvAUtils::Get_pn_reco_C_MINERvAPS(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true);
break;
case (kDalphaT):
- fXVar = FitUtils::Get_STV_dalphat_protonps(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) * (180.0 / M_PI);
+ fXVar = MINERvAUtils::Get_STV_dalphat_MINERvAPS(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) * (180.0 / M_PI);
break;
case (kDpT):
- fXVar = FitUtils::Get_STV_dpt_protonps(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) / 1000.0;
+ fXVar = MINERvAUtils::Get_STV_dpt_MINERvAPS(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) / 1000.0;
break;
case (kDphiT):
- fXVar = FitUtils::Get_STV_dphit_protonps(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) * (180.0 / M_PI);
+ fXVar = MINERvAUtils::Get_STV_dphit_MINERvAPS(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/MINERvA/MINERvA_SignalDef.cxx b/src/MINERvA/MINERvA_SignalDef.cxx
index 59f385f..7511dc8 100644
--- a/src/MINERvA/MINERvA_SignalDef.cxx
+++ b/src/MINERvA/MINERvA_SignalDef.cxx
@@ -1,471 +1,512 @@
// 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 .
-*******************************************************************************/
+ * This file is part of NUISANCE.
+ *
+ * NUISANCE is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * NUISANCE is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with NUISANCE. If not, see .
+ *******************************************************************************/
-#include "SignalDef.h"
#include "FitUtils.h"
+#include "SignalDef.h"
#include "MINERvA_SignalDef.h"
+#include "MINERvAUtils.h"
namespace SignalDef {
+// *********************************
+// MINERvA CC1pi+/- signal definition (2015 release)
+// Note: There is a 2016 release which is different to this (listed below), but
+// it is CCNpi+ and has a different W cut
+// Note2: The W cut is implemented in the class implementation in MINERvA/
+// rather than here so we can draw events that don't pass the W cut (W cut is
+// 1.4 GeV)
+// Could possibly be changed for slight speed increase since less events
+// would be used
+//
+// MINERvA signal is slightly different to MiniBooNE
+//
+// Exactly one negative muon
+// Exactly one charged pion (both + and -); however, there is a Michel e-
+// requirement but UNCLEAR IF UNFOLDED OR NOT (so don't know if should be
+// applied)
+// Exactly 1 charged pion exits (so + and - charge), however, has Michel
+// electron requirement, so look for + only here?
+// No restriction on neutral pions or other mesons
+// MINERvA has unfolded and not unfolded muon phase space for 2015
+//
+// Possible issues with the data:
+// 1) pi- is allowed in signal even when Michel cut included; most pi- is
+// efficiency corrected in GENIE 2) There is a T_pi < 350 MeV cut coming from
+// requiring a stopping pion; this is efficiency corrected in GENIE 3) There is
+// a 1.5 < Enu < 10.0 cut in signal definition 4) There is an angular muon cut
+// which is sometimes efficiency corrected (why we have bool isRestricted below)
+//
+// Nice things:
+// Much data given: with and without muon angle cuts and with and without shape
+// only data + covariance
+//
+bool isCC1pip_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
+ bool isRestricted) {
// *********************************
- // MINERvA CC1pi+/- signal definition (2015 release)
- // Note: There is a 2016 release which is different to this (listed below), but
- // it is CCNpi+ and has a different W cut
- // Note2: The W cut is implemented in the class implementation in MINERvA/
- // rather than here so we can draw events that don't pass the W cut (W cut is
- // 1.4 GeV)
- // Could possibly be changed for slight speed increase since less events
- // would be used
- //
- // MINERvA signal is slightly different to MiniBooNE
- //
- // Exactly one negative muon
- // Exactly one charged pion (both + and -); however, there is a Michel e-
- // requirement but UNCLEAR IF UNFOLDED OR NOT (so don't know if should be
- // applied)
- // Exactly 1 charged pion exits (so + and - charge), however, has Michel
- // electron requirement, so look for + only here?
- // No restriction on neutral pions or other mesons
- // MINERvA has unfolded and not unfolded muon phase space for 2015
- //
- // Possible issues with the data:
- // 1) pi- is allowed in signal even when Michel cut included; most pi- is efficiency corrected in GENIE
- // 2) There is a T_pi < 350 MeV cut coming from requiring a stopping pion; this is efficiency corrected in GENIE
- // 3) There is a 1.5 < Enu < 10.0 cut in signal definition
- // 4) There is an angular muon cut which is sometimes efficiency corrected (why we have bool isRestricted below)
- //
- // Nice things:
- // Much data given: with and without muon angle cuts and with and without shape
- // only data + covariance
- //
- bool isCC1pip_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
- bool isRestricted) {
- // *********************************
-
- // Signal is both pi+ and pi-
- // WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG
- // First, make sure it's CCINC
- if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
-
- // Allow pi+/pi-
- int piPDG[] = {211, -211};
- int nLeptons = event->NumFSLeptons();
- int nPion = event->NumFSParticle(piPDG);
-
- // Check that the desired pion exists and is the only meson
- if (nPion != 1) return false;
-
- // Check that there is only one final state lepton
- if (nLeptons != 1) return false;
-
- // MINERvA released another CC1pi+ xsec without muon unfolding!
- // here the muon angle is < 20 degrees (seen in MINOS)
- TLorentzVector pnu = event->GetHMISParticle(14)->fP;
- TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
-
- if (isRestricted) {
- double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
- if (th_nu_mu >= 20) return false;
- }
- // Extract Hadronic Mass
- double hadMass = FitUtils::Wrec(pnu, pmu);
+ // Signal is both pi+ and pi-
+ // WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG
+ // First, make sure it's CCINC
+ if (!isCCINC(event, 14, EnuMin, EnuMax))
+ return false;
- // Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition.
-#ifdef __GENIE_ENABLED__
- if (event->fType == kGENIE){
- EventRecord * gevent = static_cast(event->genie_event->event);
- const Interaction * interaction = gevent->Summary();
- const Kinematics & kine = interaction->Kine();
- double Ws = kine.W (true);
- // std::cout << "Ws versus WRec = " << Ws << " vs " << hadMass << " " << kine.W(false) << std::endl;
- hadMass = Ws * 1000.0;
- }
-#endif
- if (hadMass > 1400.0) return false;
+ // Allow pi+/pi-
+ int piPDG[] = {211, -211};
+ int nLeptons = event->NumFSLeptons();
+ int nPion = event->NumFSParticle(piPDG);
- return true;
- };
+ // Check that the desired pion exists and is the only meson
+ if (nPion != 1)
+ return false;
+ // Check that there is only one final state lepton
+ if (nLeptons != 1)
+ return false;
- // Updated MINERvA 2017 Signal using Wexp and no restriction on angle
- bool isCC1pip_MINERvA_2017(FitEvent *event, double EnuMin, double EnuMax){
+ // MINERvA released another CC1pi+ xsec without muon unfolding!
+ // here the muon angle is < 20 degrees (seen in MINOS)
+ TLorentzVector pnu = event->GetHMISParticle(14)->fP;
+ TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
- // Signal is both pi+ and pi-
- // WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG
- // First, make sure it's CCINC
- if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
+ if (isRestricted) {
+ double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
+ if (th_nu_mu >= 20)
+ return false;
+ }
- // Allow pi+/pi-
- int piPDG[] = {211, -211};
- int nLeptons = event->NumFSLeptons();
- int nPion = event->NumFSParticle(piPDG);
+ // Extract Hadronic Mass
+ double hadMass = FitUtils::Wrec(pnu, pmu);
- // Check that the desired pion exists and is the only meson
- if (nPion != 1) return false;
+ // Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition.
+#ifdef __GENIE_ENABLED__
+ if (event->fType == kGENIE) {
+ EventRecord *gevent = static_cast(event->genie_event->event);
+ const Interaction *interaction = gevent->Summary();
+ const Kinematics &kine = interaction->Kine();
+ double Ws = kine.W(true);
+ // std::cout << "Ws versus WRec = " << Ws << " vs " << hadMass << " " <<
+ // kine.W(false) << std::endl;
+ hadMass = Ws * 1000.0;
+ }
+#endif
+ if (hadMass > 1400.0)
+ return false;
- // Check that there is only one final state lepton
- if (nLeptons != 1) return false;
+ return true;
+};
- // Get Muon and Lepton Kinematics
- TLorentzVector pnu = event->GetHMISParticle(14)->fP;
- TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
+// Updated MINERvA 2017 Signal using Wexp and no restriction on angle
+bool isCC1pip_MINERvA_2017(FitEvent *event, double EnuMin, double EnuMax) {
- // Extract Hadronic Mass
- double hadMass = FitUtils::Wrec(pnu, pmu);
- // Cut on 2017 data is still 1.4 GeV
- if (hadMass > 1400.0) return false;
+ // Signal is both pi+ and pi-
+ // WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG
+ // First, make sure it's CCINC
+ if (!isCCINC(event, 14, EnuMin, EnuMax))
+ return false;
- return true;
- };
+ // Allow pi+/pi-
+ int piPDG[] = {211, -211};
+ int nLeptons = event->NumFSLeptons();
+ int nPion = event->NumFSParticle(piPDG);
+
+ // Check that the desired pion exists and is the only meson
+ if (nPion != 1)
+ return false;
+
+ // Check that there is only one final state lepton
+ if (nLeptons != 1)
+ return false;
+
+ // Get Muon and Lepton Kinematics
+ TLorentzVector pnu = event->GetHMISParticle(14)->fP;
+ TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
+ // Extract Hadronic Mass
+ double hadMass = FitUtils::Wrec(pnu, pmu);
+ // Cut on 2017 data is still 1.4 GeV
+ if (hadMass > 1400.0)
+ return false;
+
+ return true;
+};
+
+// *********************************
+// MINERvA CCNpi+/- signal definition from 2016 publication
+// Different to CC1pi+/- listed above; additional has W < 1.8 GeV
+//
+// For notes on strangeness of signal definition, see CC1pip_MINERvA
+//
+// One negative muon
+// At least one charged pion
+// 1.5 < Enu < 10
+// No restrictions on pi0 or other mesons or baryons
+// W_reconstructed (ignoring initial state motion) cut at 1.8 GeV
+//
+// Also writes number of pions (nPions) if studies on this want to be done...
+bool isCCNpip_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
+ bool isRestricted, bool isWtrue) {
// *********************************
- // MINERvA CCNpi+/- signal definition from 2016 publication
- // Different to CC1pi+/- listed above; additional has W < 1.8 GeV
- //
- // For notes on strangeness of signal definition, see CC1pip_MINERvA
- //
- // One negative muon
- // At least one charged pion
- // 1.5 < Enu < 10
- // No restrictions on pi0 or other mesons or baryons
- // W_reconstructed (ignoring initial state motion) cut at 1.8 GeV
- //
- // Also writes number of pions (nPions) if studies on this want to be done...
- bool isCCNpip_MINERvA(FitEvent *event, double EnuMin,
- double EnuMax, bool isRestricted, bool isWtrue) {
- // *********************************
-
- // First, make sure it's CCINC
- if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
-
- int nLeptons = event->NumFSLeptons();
-
- // Write the number of pions to the measurement class...
- // Maybe better to just do that inside the class?
- int nPions = event->NumFSParticle(PhysConst::pdg_charged_pions);
-
- // Check that there is a pion!
- if (nPions == 0) return false;
-
- // Check that there is only one final state lepton
- if (nLeptons != 1) return false;
-
- // Need the muon and the neutrino to check angles and W
-
- TLorentzVector pnu = event->GetNeutrinoIn()->fP;
- TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
- // MINERvA released some data with restricted muon angle
- // Here the muon angle is < 20 degrees (seen in MINOS)
- if (isRestricted) {
-
- double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
- if (th_nu_mu >= 20.) return false;
- }
- // Lastly check the W cut (always at 1.8 GeV)
- double Wrec = FitUtils::Wrec(pnu, pmu) + 0.;
+ // First, make sure it's CCINC
+ if (!isCCINC(event, 14, EnuMin, EnuMax))
+ return false;
- // Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition.
- if (isWtrue){
-#ifdef __GENIE_ENABLED__
- if (event->fType == kGENIE){
- GHepRecord* ghep = static_cast(event->genie_event->event);
- const Interaction * interaction = ghep->Summary();
- const Kinematics & kine = interaction->Kine();
- double Ws = kine.W (true);
- Wrec = Ws * 1000.0; // Say Wrec is Ws
- }
-#endif
- }
+ int nLeptons = event->NumFSLeptons();
- if (Wrec > 1800. || Wrec < 0.0) return false;
+ // Write the number of pions to the measurement class...
+ // Maybe better to just do that inside the class?
+ int nPions = event->NumFSParticle(PhysConst::pdg_charged_pions);
+ // Check that there is a pion!
+ if (nPions == 0)
+ return false;
- return true;
- };
+ // Check that there is only one final state lepton
+ if (nLeptons != 1)
+ return false;
- //********************************************************************
- bool isCCQEnumu_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
- bool fullphasespace) {
- //********************************************************************
+ // Need the muon and the neutrino to check angles and W
- if (!isCCQELike(event, 14, EnuMin, EnuMax)) return false;
+ TLorentzVector pnu = event->GetNeutrinoIn()->fP;
+ TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
+ // MINERvA released some data with restricted muon angle
+ // Here the muon angle is < 20 degrees (seen in MINOS)
+ if (isRestricted) {
- TLorentzVector pnu = event->GetHMISParticle(14)->fP;
- TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
+ double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
+ if (th_nu_mu >= 20.)
+ return false;
+ }
- double ThetaMu = pnu.Vect().Angle(pmu.Vect());
- double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 34., true);
+ // Lastly check the W cut (always at 1.8 GeV)
+ double Wrec = FitUtils::Wrec(pnu, pmu) + 0.;
- // If Restricted phase space
- if (!fullphasespace && ThetaMu > 0.34906585) return false;
+ // Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition.
+ if (isWtrue) {
+#ifdef __GENIE_ENABLED__
+ if (event->fType == kGENIE) {
+ GHepRecord *ghep = static_cast(event->genie_event->event);
+ const Interaction *interaction = ghep->Summary();
+ const Kinematics &kine = interaction->Kine();
+ double Ws = kine.W(true);
+ Wrec = Ws * 1000.0; // Say Wrec is Ws
+ }
+#endif
+ }
- // restrict energy range
- if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false;
+ if (Wrec > 1800. || Wrec < 0.0)
+ return false;
- return true;
- };
+ return true;
+};
+//********************************************************************
+bool isCCQEnumu_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
+ bool fullphasespace) {
//********************************************************************
- bool isCCQEnumubar_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
- bool fullphasespace) {
- //********************************************************************
- if (!isCCQELike(event, -14, EnuMin, EnuMax)) return false;
+ if (!isCCQELike(event, 14, EnuMin, EnuMax))
+ return false;
- TLorentzVector pnu = event->GetHMISParticle(-14)->fP;
- TLorentzVector pmu = event->GetHMFSParticle(-13)->fP;
+ TLorentzVector pnu = event->GetHMISParticle(14)->fP;
+ TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
- double ThetaMu = pnu.Vect().Angle(pmu.Vect());
- double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 30., true);
+ double ThetaMu = pnu.Vect().Angle(pmu.Vect());
+ double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 34., true);
- // If Restricted phase space
- if (!fullphasespace && ThetaMu > 0.34906585) return false;
+ // If Restricted phase space
+ if (!fullphasespace && ThetaMu > 0.34906585)
+ return false;
- // restrict energy range
- if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false;
+ // restrict energy range
+ if (Enu_rec < EnuMin || Enu_rec > EnuMax)
+ return false;
- return true;
- }
+ return true;
+};
+//********************************************************************
+bool isCCQEnumubar_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
+ bool fullphasespace) {
//********************************************************************
- bool isCCincLowRecoil_MINERvA(FitEvent *event, double EnuMin, double EnuMax) {
- //********************************************************************
-
- if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
- // Need at least one muon
- if (event->NumFSParticle(13) < 1) return false;
- TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
- TLorentzVector pnu = event->GetHMISParticle(14)->fP;
-
- // Cut on muon angle greated than 20deg
- if (cos(pnu.Vect().Angle(pmu.Vect())) < 0.93969262078) return false;
+ if (!isCCQELike(event, -14, EnuMin, EnuMax))
+ return false;
- // Cut on muon energy < 1.5 GeV
- if (pmu.E()/1000.0 < 1.5) return false;
+ TLorentzVector pnu = event->GetHMISParticle(-14)->fP;
+ TLorentzVector pmu = event->GetHMFSParticle(-13)->fP;
- return true;
- }
+ double ThetaMu = pnu.Vect().Angle(pmu.Vect());
+ double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 30., true);
+ // If Restricted phase space
+ if (!fullphasespace && ThetaMu > 0.34906585)
+ return false;
- // Used in 2014 muon+proton analysis
- // Events with muon angles up to 70 degrees
- // One right sign muon, at least one proton, no pions
- // proton kinetic energies greater than 100 MeV
- bool isCC0pi1p_MINERvA(FitEvent *event, double enumin, double enumax) {
- bool signal = (
- isCC0pi(event, 14, enumin, enumax) && // Require numu CC0pi event
- HasProtonKEAboveThreshold(event, 110.0) && // With proton above threshold
- (event->GetHMFSMuon())->P3().Angle(
- (event->GetNeutrinoIn())->P3())*180./M_PI < 70 // And muon within production angle
- );
+ // restrict energy range
+ if (Enu_rec < EnuMin || Enu_rec > EnuMax)
+ return false;
- return signal;
- }
+ return true;
+}
- // 2015 analysis just asks for 1pi0 and no pi+/pi-
- bool isCC1pi0_MINERvA_2015(FitEvent *event, double EnuMin, double EnuMax) {
- bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax);
- return CC1pi0_anu;
- }
+//********************************************************************
+bool isCCincLowRecoil_MINERvA(FitEvent *event, double EnuMin, double EnuMax) {
+ //********************************************************************
- // 2016 analysis just asks for 1pi0 and no other charged tracks. Half-open to interpretation: we go with "charged tracks" meaning pions. You'll be forgiven for thinking proton tracks should be included here too but we checked with MINERvA
- bool isCC1pi0_MINERvA_2016(FitEvent *event, double EnuMin, double EnuMax) {
- bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax);
+ if (!isCCINC(event, 14, EnuMin, EnuMax))
+ return false;
- /*
- // Additionally look for charged proton track
- // Comment: This is _NOT_ in the signal definition but was tested
- bool HasProton = event->HasFSParticle(2212);
+ // Need at least one muon
+ if (event->NumFSParticle(13) < 1)
+ return false;
+ TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
+ TLorentzVector pnu = event->GetHMISParticle(14)->fP;
- if (CC1pi0_anu) {
- if (!HasProton) {
- return true;
- } else {
+ // Cut on muon angle greated than 20deg
+ if (cos(pnu.Vect().Angle(pmu.Vect())) < 0.93969262078)
return false;
- }
- } else {
+
+ // Cut on muon energy < 1.5 GeV
+ if (pmu.E() / 1000.0 < 1.5)
return false;
- }
- */
- return CC1pi0_anu;
+ return true;
+}
+
+// Used in 2014 muon+proton analysis
+// Events with muon angles up to 70 degrees
+// One right sign muon, at least one proton, no pions
+// proton kinetic energies greater than 100 MeV
+bool isCC0pi1p_MINERvA(FitEvent *event, double enumin, double enumax) {
+ bool signal =
+ (isCC0pi(event, 14, enumin, enumax) && // Require numu CC0pi event
+ HasProtonKEAboveThreshold(event, 110.0) && // With proton above threshold
+ (event->GetHMFSMuon())->P3().Angle((event->GetNeutrinoIn())->P3()) *
+ 180. / M_PI <
+ 70 // And muon within production angle
+ );
+
+ return signal;
+}
+
+// 2015 analysis just asks for 1pi0 and no pi+/pi-
+bool isCC1pi0_MINERvA_2015(FitEvent *event, double EnuMin, double EnuMax) {
+ bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax);
+ return CC1pi0_anu;
+}
+
+// 2016 analysis just asks for 1pi0 and no other charged tracks. Half-open to
+// interpretation: we go with "charged tracks" meaning pions. You'll be forgiven
+// for thinking proton tracks should be included here too but we checked with
+// MINERvA
+bool isCC1pi0_MINERvA_2016(FitEvent *event, double EnuMin, double EnuMax) {
+ bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax);
+
+ /*
+ // Additionally look for charged proton track
+ // Comment: This is _NOT_ in the signal definition but was tested
+ bool HasProton = event->HasFSParticle(2212);
+
+ if (CC1pi0_anu) {
+ if (!HasProton) {
+ return true;
+ } else {
+ return false;
}
-
- // 2016 analysis just asks for 1pi0 and no other charged tracks
- bool isCC1pi0_MINERvA_nu(FitEvent *event, double EnuMin, double EnuMax) {
- bool CC1pi0_nu = SignalDef::isCC1pi(event, 14, 111, EnuMin, EnuMax);
- return CC1pi0_nu;
+ } else {
+ return false;
}
+ */
+
+ return CC1pi0_anu;
+}
+
+// 2016 analysis just asks for 1pi0 and no other charged tracks
+bool isCC1pi0_MINERvA_nu(FitEvent *event, double EnuMin, double EnuMax) {
+ bool CC1pi0_nu = SignalDef::isCC1pi(event, 14, 111, EnuMin, EnuMax);
+ return CC1pi0_nu;
+}
+//********************************************************************
+bool isCC0pi_MINERvAPTPZ(FitEvent *event, int nuPDG, double emin, double emax) {
//********************************************************************
- bool isCC0pi_MINERvAPTPZ(FitEvent* event, int nuPDG, double emin, double emax){
- //********************************************************************
- // Check it's CCINC
- if (!SignalDef::isCCINC(event, nuPDG, emin, emax)) return false;
-
- // Make Angle Cut > 20.0
- TLorentzVector pnu = event->GetHMISParticle(14)->fP;
- TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
- double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
- if (th_nu_mu >= 20.0) return false;
-
- int genie_n_muons = 0;
- int genie_n_mesons = 0;
- int genie_n_heavy_baryons_plus_pi0s = 0;
- int genie_n_photons = 0;
-
- for(unsigned int i = 0; i < event->NParticles(); ++i) {
- FitParticle* p = event->GetParticle(i);
- if (p->Status() != kFinalState) continue;
-
- int pdg = p->fPID;
- double energy = p->fP.E();
-
- if( pdg == 13 ) {
- genie_n_muons++;
- } else if( pdg == 22 && energy > 10.0 ) {
- genie_n_photons++;
- } else if( abs(pdg) == 211 || abs(pdg) == 321 || abs(pdg) == 323 || pdg == 111 || pdg == 130 || pdg == 310 || pdg == 311 || pdg == 313 || abs(pdg) == 221 || abs(pdg) == 331 ) {
- genie_n_mesons++;
- } else if( pdg == 3112 || pdg == 3122 || pdg == 3212 || pdg == 3222 || pdg == 4112 ||
- pdg == 4122 || pdg == 4212 || pdg == 4222 || pdg == 411 || pdg == 421 || pdg == 111 ) {
- genie_n_heavy_baryons_plus_pi0s++;
- }
- }
+ // Check it's CCINC
+ if (!SignalDef::isCCINC(event, nuPDG, emin, emax))
+ return false;
- if( genie_n_muons == 1 &&
- genie_n_mesons == 0 &&
- genie_n_heavy_baryons_plus_pi0s == 0 &&
- genie_n_photons == 0 ) return true;
+ // Make Angle Cut > 20.0
+ TLorentzVector pnu = event->GetHMISParticle(14)->fP;
+ TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
+ double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
+ if (th_nu_mu >= 20.0)
return false;
+
+ int genie_n_muons = 0;
+ int genie_n_mesons = 0;
+ int genie_n_heavy_baryons_plus_pi0s = 0;
+ int genie_n_photons = 0;
+
+ for (unsigned int i = 0; i < event->NParticles(); ++i) {
+ FitParticle *p = event->GetParticle(i);
+ if (p->Status() != kFinalState)
+ continue;
+
+ int pdg = p->fPID;
+ double energy = p->fP.E();
+
+ if (pdg == 13) {
+ genie_n_muons++;
+ } else if (pdg == 22 && energy > 10.0) {
+ genie_n_photons++;
+ } else if (abs(pdg) == 211 || abs(pdg) == 321 || abs(pdg) == 323 ||
+ pdg == 111 || pdg == 130 || pdg == 310 || pdg == 311 ||
+ pdg == 313 || abs(pdg) == 221 || abs(pdg) == 331) {
+ genie_n_mesons++;
+ } else if (pdg == 3112 || pdg == 3122 || pdg == 3212 || pdg == 3222 ||
+ pdg == 4112 || pdg == 4122 || pdg == 4212 || pdg == 4222 ||
+ pdg == 411 || pdg == 421 || pdg == 111) {
+ genie_n_heavy_baryons_plus_pi0s++;
+ }
}
- // **************************************************
- // Section VI Event Selection of https://journals.aps.org/prd/pdf/10.1103/PhysRevD.97.052002
- // Anti-neutrino charged-current
- // Post-FSI final states without
- // mesons,
- // prompt photons above nuclear deexcitation energies
- // heavy baryons
- // protons above kinetic energy of 120 MeV
- // Muon-neutrino angle of 20 degrees
- // Parallel muon momentum: 1.5 GeV < P|| < 15 GeV --- N.B. APPARENTLY NOT INCLUDED, see below
- // Transverse muon momentum: pT < 1.5 GeV --- N.B. APPARENTLY NOT INCLUDED, see below
- bool isCC0pi_anti_MINERvAPTPZ(FitEvent* event, int nuPDG, double emin, double emax){
+ if (genie_n_muons == 1 && genie_n_mesons == 0 &&
+ genie_n_heavy_baryons_plus_pi0s == 0 && genie_n_photons == 0)
+ return true;
+ return false;
+}
+
+// **************************************************
+// Section VI Event Selection of
+// https://journals.aps.org/prd/pdf/10.1103/PhysRevD.97.052002 Anti-neutrino
+// charged-current Post-FSI final states without
+// mesons,
+// prompt photons above nuclear deexcitation
+// energies heavy baryons protons above kinetic
+// energy of 120 MeV
+// Muon-neutrino angle of 20 degrees
+// Parallel muon momentum: 1.5 GeV < P|| < 15 GeV --- N.B. APPARENTLY NOT
+// INCLUDED, see below Transverse muon momentum: pT < 1.5 GeV --- N.B.
+// APPARENTLY NOT INCLUDED, see below
+bool isCC0pi_anti_MINERvAPTPZ(FitEvent *event, int nuPDG, double emin,
+ double emax) {
// **************************************************
- // Check it's CCINC
- if (!SignalDef::isCCINC(event, nuPDG, emin, emax)) return false;
- TLorentzVector pnu = event->GetNeutrinoIn()->fP;
- TLorentzVector pmu = event->GetHMFSParticle(-13)->fP;
- // Make Angle Cut > 20.0
- double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
- if (th_nu_mu >= 20.0) return false;
-
- // Heidi Schellman (schellmh@science.oregonstate.edu) assured me that the p_t and p_z (or p_||) cuts aren't actually implemented as a signal definition: they're only implemented in the binning for p_t and p_z (but not Q2QE and EnuQE)
- /*
- // Cut on pT and pZ
- Double_t px = pmu.X()/1.E3;
- Double_t py = pmu.Y()/1.E3;
- Double_t pt = sqrt(px*px+py*py);
-
- // Don't want to assume the event generators all have neutrino coming along z
- // pz is muon momentum projected onto the neutrino direction
- Double_t pz = pmu.Vect().Dot(pnu.Vect()*(1.0/pnu.Vect().Mag()))/1.E3;
- if (pz > 15 || pz < 1.5) return false;
- if (pt > 1.5) return false;
- */
-
- // Find if there are any protons above 120 MeV kinetic energy
- if (HasProtonKEAboveThreshold(event, 120.0)) return false;
-
- // Particle counters
- int genie_n_muons = 0;
- int genie_n_mesons = 0;
- int genie_n_heavy_baryons_plus_pi0s = 0;
- int genie_n_photons = 0;
- // Loop over the particles in the event and count them up
- for(unsigned int i = 0; i < event->NParticles(); ++i) {
- FitParticle* p = event->GetParticle(i);
- if (p->Status() != kFinalState) continue;
-
- int pdg = p->fPID;
- double energy = p->fP.E();
-
- // Any charged muons
- if( abs(pdg) == 13 ) {
- genie_n_muons++;
- // De-excitation photons
- } else if( pdg == 22 && energy > 10.0 ) {
- genie_n_photons++;
- // Mesons
- } else if( abs(pdg) == 211 || abs(pdg) == 321 || abs(pdg) == 323 || pdg == 111 || pdg == 130 || pdg == 310 || pdg == 311 || pdg == 313 || abs(pdg) == 221 || abs(pdg) == 331 ) {
- genie_n_mesons++;
- // Heavy baryons and pi0s
- } else if( abs(pdg) == 3112 || abs(pdg) == 3122 || abs(pdg) == 3212 || abs(pdg) == 3222 || abs(pdg) == 4112 ||
- abs(pdg) == 4122 || abs(pdg) == 4212 || abs(pdg) == 4222 || abs(pdg) == 411 || abs(pdg) == 421 ||
- abs(pdg) == 111 ) {
- genie_n_heavy_baryons_plus_pi0s++;
- }
- }
+ // Check it's CCINC
+ if (!SignalDef::isCCINC(event, nuPDG, emin, emax))
+ return false;
+ TLorentzVector pnu = event->GetNeutrinoIn()->fP;
+ TLorentzVector pmu = event->GetHMFSParticle(-13)->fP;
+ // Make Angle Cut > 20.0
+ double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
+ if (th_nu_mu >= 20.0)
+ return false;
- // Look for one muon with no mesons, heavy baryons or deexcitation photons
- if( genie_n_muons == 1 && genie_n_mesons == 0 && genie_n_heavy_baryons_plus_pi0s == 0 && genie_n_photons == 0 ) return true;
+ // Heidi Schellman (schellmh@science.oregonstate.edu) assured me that the p_t
+ // and p_z (or p_||) cuts aren't actually implemented as a signal definition:
+ // they're only implemented in the binning for p_t and p_z (but not Q2QE and
+ // EnuQE)
+ /*
+ // Cut on pT and pZ
+ Double_t px = pmu.X()/1.E3;
+ Double_t py = pmu.Y()/1.E3;
+ Double_t pt = sqrt(px*px+py*py);
+
+ // Don't want to assume the event generators all have neutrino coming along z
+ // pz is muon momentum projected onto the neutrino direction
+ Double_t pz = pmu.Vect().Dot(pnu.Vect()*(1.0/pnu.Vect().Mag()))/1.E3;
+ if (pz > 15 || pz < 1.5) return false;
+ if (pt > 1.5) return false;
+ */
+
+ // Find if there are any protons above 120 MeV kinetic energy
+ if (HasProtonKEAboveThreshold(event, 120.0))
return false;
- }
- // MINERvA CC0pi transverse variables signal defintion
- bool isCC0piNp_MINERvA_STV(FitEvent *event, double EnuMin, double EnuMax) {
-
- // Require a numu CC0pi event
- if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false;
-
- // Require at least one FS proton
- if (event->NumFSParticle(2212) == 0) return false;
-
- TLorentzVector pnu = event->GetHMISParticle(14)->fP;
- TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
-
- // Muon momentum cuts
- if (pmu.Vect().Mag() < 1500 || pmu.Vect().Mag() > 10000) return false;
- // Muon angle cuts
- if (pmu.Vect().Angle(pnu.Vect()) > (M_PI/180.0)*20.0) return false;
-
- // Since there are emany protons allowed we can't just use the highest proton momentum candidate and place cuts on it
- // Get the stack of protons
- std::vector Protons = event->GetAllFSProton();
- // Count how many protons pass the threshold
- int nProtonsAboveThreshold = 0;
- for (size_t i = 0; i < Protons.size(); ++i) {
- if (Protons[i]->p() > 450 &&
- Protons[i]->p() < 1200 &&
- Protons[i]->P3().Angle(pnu.Vect()) < (M_PI/180.0)*70.0) {
- nProtonsAboveThreshold++;
- }
+ // Particle counters
+ int genie_n_muons = 0;
+ int genie_n_mesons = 0;
+ int genie_n_heavy_baryons_plus_pi0s = 0;
+ int genie_n_photons = 0;
+ // Loop over the particles in the event and count them up
+ for (unsigned int i = 0; i < event->NParticles(); ++i) {
+ FitParticle *p = event->GetParticle(i);
+ if (p->Status() != kFinalState)
+ continue;
+
+ int pdg = p->fPID;
+ double energy = p->fP.E();
+
+ // Any charged muons
+ if (abs(pdg) == 13) {
+ genie_n_muons++;
+ // De-excitation photons
+ } else if (pdg == 22 && energy > 10.0) {
+ genie_n_photons++;
+ // Mesons
+ } else if (abs(pdg) == 211 || abs(pdg) == 321 || abs(pdg) == 323 ||
+ pdg == 111 || pdg == 130 || pdg == 310 || pdg == 311 ||
+ pdg == 313 || abs(pdg) == 221 || abs(pdg) == 331) {
+ genie_n_mesons++;
+ // Heavy baryons and pi0s
+ } else if (abs(pdg) == 3112 || abs(pdg) == 3122 || abs(pdg) == 3212 ||
+ abs(pdg) == 3222 || abs(pdg) == 4112 || abs(pdg) == 4122 ||
+ abs(pdg) == 4212 || abs(pdg) == 4222 || abs(pdg) == 411 ||
+ abs(pdg) == 421 || abs(pdg) == 111) {
+ genie_n_heavy_baryons_plus_pi0s++;
}
+ }
- if (nProtonsAboveThreshold == 0) return false;
+ // Look for one muon with no mesons, heavy baryons or deexcitation photons
+ if (genie_n_muons == 1 && genie_n_mesons == 0 &&
+ genie_n_heavy_baryons_plus_pi0s == 0 && genie_n_photons == 0)
return true;
- };
+ return false;
+}
+
+// MINERvA CC0pi transverse variables signal defintion
+bool isCC0piNp_MINERvA_STV(FitEvent *event, double EnuMin, double EnuMax) {
+
+ // Require a numu CC0pi event
+ if (!isCC0pi(event, 14, EnuMin, EnuMax))
+ return false;
+
+ // Require at least one FS proton
+ if (event->NumFSParticle(2212) == 0)
+ return false;
+
+ TLorentzVector pnu = event->GetHMISParticle(14)->fP;
+ TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
+
+ // Muon momentum cuts
+ if (pmu.Vect().Mag() < 1500 || pmu.Vect().Mag() > 10000)
+ return false;
+ // Muon angle cuts
+ if (pmu.Vect().Angle(pnu.Vect()) > (M_PI / 180.0) * 20.0)
+ return false;
+
+ const double ctheta_cut = cos((M_PI / 180.0) * 70.0);
+ //Did you find a proton in the PS?
+ if (MINERvAUtils::GetProtonInRange(event, 450, 1200, ctheta_cut).E() == 0)
+ return false;
+ return true;
+};
-} // end namespace
+} // namespace SignalDef
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
index d322a90..05bd3cc 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
@@ -1,94 +1,90 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "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_protonps(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true);
+ fXVar = FitUtils::Get_STV_dalphat_HMProton(event, 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_1Ddat_nu.h b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.h
index f226217..2387936 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.h
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.h
@@ -1,40 +1,35 @@
// 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 .
*******************************************************************************/
#ifndef T2K_CC0pinp_STV_XSec_1Ddat_nu_H_SEEN
#define T2K_CC0pinp_STV_XSec_1Ddat_nu_H_SEEN
#include "Measurement1D.h"
class T2K_CC0pinp_STV_XSec_1Ddat_nu : public Measurement1D {
public:
T2K_CC0pinp_STV_XSec_1Ddat_nu(nuiskey samplekey);
virtual ~T2K_CC0pinp_STV_XSec_1Ddat_nu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
- private:
- // Proton cuts
- double ProtonMinCut;
- double ProtonMaxCut;
- double ProtonCosThetaCut;
};
#endif
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
index bded9c4..a01ed78 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
@@ -1,93 +1,90 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "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_protonps(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true);
+ fXVar = FitUtils::Get_STV_dphit_HMProton(event, 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_1Ddphit_nu.h b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.h
index 8f21503..de4ab66 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.h
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.h
@@ -1,40 +1,35 @@
// 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 .
*******************************************************************************/
#ifndef T2K_CC0pinp_STV_XSec_1Ddphit_nu_H_SEEN
#define T2K_CC0pinp_STV_XSec_1Ddphit_nu_H_SEEN
#include "Measurement1D.h"
class T2K_CC0pinp_STV_XSec_1Ddphit_nu : public Measurement1D {
public:
T2K_CC0pinp_STV_XSec_1Ddphit_nu(nuiskey samplekey);
virtual ~T2K_CC0pinp_STV_XSec_1Ddphit_nu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
- private:
- // Proton cuts
- double ProtonMinCut;
- double ProtonMaxCut;
- double ProtonCosThetaCut;
};
#endif
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
index b44be72..e61d59b 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
@@ -1,93 +1,90 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "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_protonps(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true) / 1000.0;
+ fXVar = FitUtils::Get_STV_dpt_HMProton(event, 14, true)*1E-3;
return;
};
//********************************************************************
bool T2K_CC0pinp_STV_XSec_1Ddpt_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.h b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.h
index d828684..a83e260 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.h
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.h
@@ -1,40 +1,35 @@
// 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 .
*******************************************************************************/
#ifndef T2K_CC0pinp_STV_XSec_1Ddpt_nu_H_SEEN
#define T2K_CC0pinp_STV_XSec_1Ddpt_nu_H_SEEN
#include "Measurement1D.h"
class T2K_CC0pinp_STV_XSec_1Ddpt_nu : public Measurement1D {
public:
T2K_CC0pinp_STV_XSec_1Ddpt_nu(nuiskey samplekey);
virtual ~T2K_CC0pinp_STV_XSec_1Ddpt_nu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
- private:
- // Proton cuts
- double ProtonMinCut;
- double ProtonMaxCut;
- double ProtonCosThetaCut;
};
#endif
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx
index 706afb9..19e493d 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx
@@ -1,110 +1,111 @@
#include
#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) {
// Sample overview ---------------------------------------------------
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");
// 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");
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;
// Reconstruct the neutrino
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
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);
// 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
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()));
// 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());
// And the x-axis is then simply perpendicular to z and x
TVector3 xAxis = yAxis.Cross(zAxis);
xAxis *= 1.0 / double(xAxis.Mag());
double x = Ppip.Vect().Dot(xAxis);
double y = Ppip.Vect().Dot(yAxis);
double newphi = atan2(y, x);
fXVar = newphi;
return;
};
//********************************************************************
bool T2K_CC1pip_CH_XSec_1DAdlerPhi_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCC1pip_T2K_arxiv1909_03936(
- event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionHighEff);
+ event, EnuMin, EnuMax,
+ SignalDef::kMuonHighEff | SignalDef::kPionVFwd | SignalDef::kPionHighMom);
}
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx
index db554db..8b1ccac 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx
@@ -1,101 +1,102 @@
#include
#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) {
// Sample overview ---------------------------------------------------
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");
// 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");
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;
// 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 Ppip = event->GetHMFSParticle(211)->fP;
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);
// 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
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()));
// 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_arxiv1909_03936(
- event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionHighEff);
+ event, EnuMin, EnuMax,
+ SignalDef::kMuonHighEff | SignalDef::kPionVFwd | SignalDef::kPionHighMom);
}
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx
index 0c16689..5da96bc 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx
@@ -1,99 +1,102 @@
#include
#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_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");
+ 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");
-
+ 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)
+ 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 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;
+ // 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) {
-//********************************************************************
+ //********************************************************************
return SignalDef::isCC1pip_T2K_arxiv1909_03936(
- event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionHighEff);
+ event, EnuMin, EnuMax,
+ SignalDef::kMuonHighEff | SignalDef::kPionVFwd | SignalDef::kPionHighMom);
}
-
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx
index 64a5dc6..9c7c211 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx
@@ -1,70 +1,71 @@
#include
-#include "T2K_SignalDef.h"
#include "T2K_CC1pip_CH_XSec_1Dppi_nu.h"
-
+#include "T2K_SignalDef.h"
//********************************************************************
T2K_CC1pip_CH_XSec_1Dppi_nu::T2K_CC1pip_CH_XSec_1Dppi_nu(nuiskey samplekey) {
-//********************************************************************
+ //********************************************************************
// Sample overview ---------------------------------------------------
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");
+ 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");
-
+ 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)
+ 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 ppip = FitUtils::p(Ppip);
fXVar = ppip;
return;
};
//********************************************************************
bool T2K_CC1pip_CH_XSec_1Dppi_nu::isSignal(FitEvent *event) {
-//********************************************************************
+ //********************************************************************
return SignalDef::isCC1pip_T2K_arxiv1909_03936(
- event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionFwdHighMom);
+ event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionVFwd);
}
-
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx
index fe6c19a..ed4002c 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx
@@ -1,71 +1,72 @@
#include
#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) {
// Sample overview ---------------------------------------------------
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");
// 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");
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)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double thmupi = FitUtils::th(Pmu, Ppip);
fXVar = thmupi;
// std::cout << thmupi << std::endl;
return;
};
//********************************************************************
bool T2K_CC1pip_CH_XSec_1Dthmupi_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCC1pip_T2K_arxiv1909_03936(
- event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionHighEff);
+ event, EnuMin, EnuMax,
+ SignalDef::kMuonHighEff | SignalDef::kPionVFwd | SignalDef::kPionHighMom);
}
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx
index 27b25cb..9028886 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx
@@ -1,66 +1,67 @@
#include
#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_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");
// 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");
SetShapeCovar();
fDataHist->Scale(1E-38);
FinaliseMeasurement();
};
void T2K_CC1pip_CH_XSec_1Dthpi_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(211) == 0)
return;
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) {
//********************************************************************
return SignalDef::isCC1pip_T2K_arxiv1909_03936(
- event, EnuMin, EnuMax, SignalDef::kMuonHighEff | SignalDef::kPionFwdHighMom);
+ event, EnuMin, EnuMax,
+ SignalDef::kMuonHighEff | SignalDef::kPionFwd | SignalDef::kPionHighMom);
}
diff --git a/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx
index b5b7c7e..161a179 100644
--- a/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx
+++ b/src/T2K/T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx
@@ -1,213 +1,215 @@
#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) {
//********************************************************************
// Sample overview ---------------------------------------------------
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.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.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) /
double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
// SetDataValues( fSettings.GetDataInput() );
// SetCovarMatrix( fSettings.GetCovarInput() );
SetHistograms();
- //fFullCovar = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(),
- //"Covariance_pmu_thetamu");
+ // fFullCovar = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(),
+ //"Covariance_pmu_thetamu");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
SetShapeCovar();
/*
for (int i = 0; i < covar->GetNrows(); ++i) {
for (int j = 0; j < covar->GetNrows(); ++j) {
- if (i == j) std::cout << i << " " << j << " = " << 1/sqrt((*covar)(i,j)) << std::endl;
+ if (i == j) std::cout << i << " " << j << " = " << 1/sqrt((*covar)(i,j))
+ << std::endl;
}
}
throw;
*/
// 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 = "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());
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);
SetAutoProcessTH1(fMCHist_Slices[i]);
fMCHist_Slices[i]->Reset();
fMCHist_Slices[i]->SetLineColor(kRed);
- //nbins += slice->GetXaxis()->GetNbins();
- nbins += slice->GetXaxis()->GetNbins()-1;
+ // nbins += slice->GetXaxis()->GetNbins();
+ nbins += slice->GetXaxis()->GetNbins() - 1;
}
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()-1; ++j) {
+ for (int j = 0; j < fDataHist_Slices[i]->GetXaxis()->GetNbins() - 1; ++j) {
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));
fDataHist->GetXaxis()->SetBinLabel(bincount, title);
bincount++;
}
}
fDataHist->GetXaxis()->SetTitle(fSettings.GetS("xtitle").c_str());
fDataHist->GetYaxis()->SetTitle(fSettings.GetS("ytitle").c_str());
// Get the covariance
- TMatrixDSym *temp = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(), "Covariance_pmu_thetamu");
+ TMatrixDSym *temp = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(),
+ "Covariance_pmu_thetamu");
int ncovbins = temp->GetNrows();
- fFullCovar = new TMatrixDSym(ncovbins-4);
+ fFullCovar = new TMatrixDSym(ncovbins - 4);
if (ncovbins != fDataHist->GetXaxis()->GetNbins()) {
NUIS_ERR(FTL, "Number of bins in covariance matrix does not match data");
}
// Number of costhetamu slices is nslices
- // Number of pmu slices is
+ // Number of pmu slices is
int count1 = 0;
- for (int i = 0; i < ncovbins-4; ++i) {
+ for (int i = 0; i < ncovbins - 4; ++i) {
int count2 = 0;
- for (int j = 0; j < ncovbins-4; ++j) {
+ for (int j = 0; j < ncovbins - 4; ++j) {
// 1E79 matched to diagonal error
(*fFullCovar)(count1, count2) = (*temp)(i, j);
count2++;
}
count1++;
}
// Now reorganise the rows
delete temp;
};
void T2K_CC1pip_CH_XSec_2Dpmucosmu_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(13) == 0)
return;
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() {
const int nslices = 4;
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()-1; j++) {
+ for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX() - 1; 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) {
// 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) {
//********************************************************************
return SignalDef::isCC1pip_T2K_arxiv1909_03936(event, EnuMin, EnuMax,
- SignalDef::kMuonFwd);
+ SignalDef::kMuonFwd);
};
diff --git a/src/T2K/T2K_SignalDef.cxx b/src/T2K/T2K_SignalDef.cxx
index 24edfa1..62ebf01 100644
--- a/src/T2K/T2K_SignalDef.cxx
+++ b/src/T2K/T2K_SignalDef.cxx
@@ -1,342 +1,331 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "T2K_SignalDef.h"
#include "FitUtils.h"
namespace SignalDef {
// T2K H2O signal definition
// 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;
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
//
// 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;
}
TLorentzVector Pnu = event->GetHMISParticle(14)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double cos_th_mu = cos(FitUtils::th(Pnu, Pmu));
if (pscuts == kMuonFwd) {
return (cos_th_mu > 0);
}
double p_mu = FitUtils::p(Pmu) * 1000;
if (pscuts & kMuonHighEff) {
if ((cos_th_mu <= 0.2) || (p_mu <= 200)) {
return false;
}
}
- 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 & kPionFwd) && (cos_th_pi <= 0)) {
+ return false;
}
- if (pscuts & kPionHighMom) {
- return (p_pi > 200);
+ if ((pscuts & kPionVFwd) && (cos_th_pi <= 0.2)) {
+ return false;
}
- if (pscuts & kMuonHighEff) {
- return ((cos_th_pi > 0.0) && (p_pi > 200));
+ if ((pscuts & kPionHighMom) && (p_pi <= 200)) {
+ return false;
}
- return false;
+ return true;
};
bool isT2K_CC0pi(FitEvent *event, double EnuMin, double EnuMax, int ana) {
// Require a numu CC0pi event
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 (ana == kAnalysis_II) {
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;
// Require at least one FS proton
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;
// Proton phase space
if (pp.Vect().Mag() < 500) {
return false;
}
// Need exactly one proton with 500 MeV or more momentum
std::vector 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;
return true;
}
// 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;
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;
return true;
}
bool isT2K_CC0piNp(FitEvent *event, double EnuMin, double EnuMax) {
// Require a numu CC0pi event
if (!isCC0pi(event, 14, EnuMin, EnuMax))
return false;
// Require at least one FS proton
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;
// Proton phase space
if (pp.Vect().Mag() < 500) {
return false;
}
// Need exactly one proton with 500 MeV or more momentum
std::vector 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;
return true;
}
bool isT2K_CC0pi0p(FitEvent *event, double EnuMin, double EnuMax) {
// Require a numu CC0pi event
if (!isCC0pi(event, 14, EnuMin, EnuMax))
return false;
// Require at least one FS proton
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;
// 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;
// Require at least one FS proton
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;
// 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;
// Require at least one FS proton
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;
// 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;
// Require at least one FS proton
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;
// Muon phase space
// 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 aec9b64..cbc90b2 100644
--- a/src/T2K/T2K_SignalDef.h
+++ b/src/T2K/T2K_SignalDef.h
@@ -1,56 +1,56 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef T2K_SIGNALDEF_H_SEEN
#define T2K_SIGNALDEF_H_SEEN
#include "SignalDef.h"
namespace SignalDef {
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)
+ kMuonFwd = (1 << 0), // cos(th_mu) > 0
+ kMuonHighEff = (1 << 1), // cos(th_mu) > 0.2, pmu > 200
+ kPionFwd = (1 << 2), // cos(th_pi) > 0
+ kPionVFwd = (1 << 3), // cos(th_pi) > 0.2
+ kPionHighMom = (1 << 4) // ppi > 200
};
bool isCC1pip_T2K_arxiv1909_03936(FitEvent *event, double EnuMin, double EnuMax,
int cuts);
enum PRD93112012_Ana {
kAnalysis_I,
kAnalysis_II,
};
bool isT2K_CC0pi(FitEvent *event, double EnuMin, double EnuMax,
int analysis);
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
diff --git a/src/Utils/FitUtils.cxx b/src/Utils/FitUtils.cxx
index 5e9d661..3a6b424 100644
--- a/src/Utils/FitUtils.cxx
+++ b/src/Utils/FitUtils.cxx
@@ -1,1158 +1,1102 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "FitUtils.h"
/*
MISC Functions
*/
//********************************************************************
double *FitUtils::GetArrayFromMap(std::vector invals,
std::map 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 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
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 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) {
return EnuQErec(pmu, cos(pnu.Vect().Angle(pmu.Vect())), binding, 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 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);
return q2qe;
}
double FitUtils::EnuQErec(double pl, double costh, double binding,
bool neutrino) {
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
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
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 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
// 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
// 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 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 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
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 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 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;
// 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;
// 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;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == 13)
continue;
// Skip Neutrons particles
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
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();
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;
int PID = event->PartInfo(i)->fPID;
// Neutrons
if (PID == 2112) {
// Adding kinetic energy of neutron
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;
}
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) {
+Double_t FitUtils::GetDeltaPhiT(TVector3 const &V_lepton,
+ TVector3 const &V_other, TVector3 const &Normal,
+ bool PiMinus) {
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 FitUtils::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) {
+Double_t FitUtils::GetDeltaAlphaT(TVector3 const &V_lepton,
+ TVector3 const &V_other,
+ TVector3 const &Normal, bool PiMinus) {
TVector3 DeltaPT = GetDeltaPT(V_lepton, V_other, Normal);
return GetDeltaPhiT(V_lepton, DeltaPT, Normal, PiMinus);
}
-double FitUtils::Get_STV_dpt_protonps(FitEvent *event, double ProtonMinCut,
- double ProtonMaxCut,
- double ProtonCosThetaCut, int ISPDG,
- bool Is0pi) {
+double FitUtils::Get_STV_dpt_HMProton(FitEvent *event, 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_MeV and
+ // ProtonMaxCut_MeV MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = event->GetHMFSParticle(2212)->fP;
// 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();
+ TLorentzVector ppi = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
+ HadronP += ppi.Vect();
}
return GetDeltaPT(LeptonP, HadronP, NuP).Mag();
}
-double FitUtils::Get_STV_dphit_protonps(FitEvent *event, double ProtonMinCut,
- double ProtonMaxCut,
- double ProtonCosThetaCut, int ISPDG,
+double FitUtils::Get_STV_dphit_HMProton(FitEvent *event, 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_MeV and
+ // ProtonMaxCut_MeV MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = event->GetHMFSParticle(2212)->fP;
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();
+ TLorentzVector ppi = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
+ HadronP += ppi.Vect();
}
return GetDeltaPhiT(LeptonP, HadronP, NuP);
}
-double FitUtils::Get_STV_dalphat_protonps(FitEvent *event, double ProtonMinCut,
- double ProtonMaxCut,
- double ProtonCosThetaCut, int ISPDG,
+double FitUtils::Get_STV_dalphat_HMProton(FitEvent *event, 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_MeV and
+ // ProtonMaxCut_MeV MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = event->GetHMFSParticle(2212)->fP;
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();
+ TLorentzVector ppi = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
+ HadronP += ppi.Vect();
}
return GetDeltaAlphaT(LeptonP, HadronP, NuP);
}
// As defined in PhysRevC.95.065501
// Using prescription from arXiv 1805.05486
// Returns in GeV
-double FitUtils::Get_pn_reco_C_protonps(FitEvent *event, double ProtonMinCut,
- double ProtonMaxCut,
- double ProtonCosThetaCut, int ISPDG,
+double FitUtils::Get_pn_reco_C_HMProton(FitEvent *event, int ISPDG,
bool Is0pi) {
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);
+ TLorentzVector Pprot = event->GetHMFSParticle(2212)->fP;
TVector3 HadronP = Pprot.Vect();
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();
+ TLorentzVector ppi = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
+ HadronP += ppi.Vect();
}
TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP);
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 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 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;
return pn_reco;
}
-double FitUtils::Get_pn_reco_Ar_protonps(FitEvent *event, double ProtonMinCut,
- double ProtonMaxCut,
- double ProtonCosThetaCut, int ISPDG,
+double FitUtils::Get_pn_reco_Ar_HMProton(FitEvent *event, int ISPDG,
bool Is0pi) {
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);
+ TLorentzVector Pprot = event->GetHMFSParticle(2212)->fP;
TVector3 HadronP = Pprot.Vect();
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();
+ TLorentzVector ppi = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
+ HadronP += ppi.Vect();
}
TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP);
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 map = ma - mn + 0.02713; // remnant mass
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 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;
-
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) {
- // Get the neutrino
- TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
- int HMFSProton = 0;
- double HighestMomentum = 0.0;
- // Get the stack of protons
- std::vector Protons = event->GetAllFSProton();
- for (size_t i = 0; i < Protons.size(); ++i) {
- 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) {
// 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
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()));
// 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) {
// 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
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()));
// 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());
// And the x-axis is then simply perpendicular to z and x
TVector3 xAxis = yAxis.Cross(zAxis);
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);
// Convert negative angles to positive
if (newphi < 0.0)
newphi += 360.0;
return newphi;
}
//********************************************************************
double FitUtils::ppInfK(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
double el = pmu.E() / 1000.;
// 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) {
//********************************************************************
// 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
double pl_x = pmu.X() / 1000.;
double pl_y = pmu.Y() / 1000.;
double pl_z = pmu.Z() / 1000.;
double enu = EnuQErec(pmu, costh, binding, neutrino);
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
double el = pmu.E() / 1000.;
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);
return cth_inf;
};
diff --git a/src/Utils/FitUtils.h b/src/Utils/FitUtils.h
index 1ba9e85..a993314 100644
--- a/src/Utils/FitUtils.h
+++ b/src/Utils/FitUtils.h
@@ -1,256 +1,197 @@
// 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 .
*******************************************************************************/
#ifndef FITUTILS_H_SEEN
#define FITUTILS_H_SEEN
#include
#include
#include
#include
#include
#include
#include "FitEvent.h"
#include "TGraph.h"
#include "TH2Poly.h"
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#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 invals,
std::map inmap);
/// Returns kinetic energy of particle
double T(TLorentzVector part);
/// Returns momentum of particle
double p(TLorentzVector part);
double p(FitParticle *part);
/// Returns angle between particles (_NOT_ cosine!)
double th(TLorentzVector part, TLorentzVector 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 const fps);
double SumTE_PartVect(std::vector 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);
/// Function returns the reconstructed E_{nu} values
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);
//! 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);
//! Function returns the reconstructed E_{nu} values
double EnuQErec(double pl, double costh, double binding, bool neutrino = true);
+Double_t GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other,
+ TVector3 const &Normal, bool PiMinus = false);
+TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other,
+ TVector3 const &Normal);
+Double_t GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other,
+ TVector3 const &Normal, bool PiMinus = false);
+double Get_STV_dpt_HMProton(FitEvent *event, int ISPDG, bool Is0pi);
+double Get_STV_dphit_HMProton(FitEvent *event, int ISPDG, bool Is0pi);
+double Get_STV_dalphat_HMProton(FitEvent *event, int ISPDG, bool Is0pi);
+double Get_pn_reco_C_HMProton(FitEvent *event, int ISPDG, bool Is0pi);
+double Get_pn_reco_Ar_HMProton(FitEvent *event, int ISPDG, bool Is0pi);
+
/*
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_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_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_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_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);
-}
-
-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);
-}
-
-// 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);
-
// 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