diff --git a/data/MINERvA/CC0pi/proton_2014/muon_Q2QE_nu_data.txt b/data/MINERvA/CC0pi/proton_2014/muon_Q2QE_nu_data.txt
new file mode 100644
index 0000000..102f353
--- /dev/null
+++ b/data/MINERvA/CC0pi/proton_2014/muon_Q2QE_nu_data.txt
@@ -0,0 +1,5 @@
+0.2 1.164
+0.4 0.534
+0.8 0.222
+1.2 0.110
+2.0 0.00
diff --git a/data/MINERvA/CCQE/proton_Q2QE_nu_covar.txt b/data/MINERvA/CC0pi/proton_2014/proton_Q2QE_nu_covar.txt
old mode 100755
new mode 100644
similarity index 100%
rename from data/MINERvA/CCQE/proton_Q2QE_nu_covar.txt
rename to data/MINERvA/CC0pi/proton_2014/proton_Q2QE_nu_covar.txt
diff --git a/data/MINERvA/CC0pi/proton_2014/proton_Q2QE_nu_covar_shape.txt b/data/MINERvA/CC0pi/proton_2014/proton_Q2QE_nu_covar_shape.txt
new file mode 100644
index 0000000..7f50e51
--- /dev/null
+++ b/data/MINERvA/CC0pi/proton_2014/proton_Q2QE_nu_covar_shape.txt
@@ -0,0 +1,7 @@
+1.00 0.52 -0.36 -0.83 -0.86 -0.63 -0.24
+0.52 1.00 0.10 -0.22 -0.33 -0.75 -0.30
+-0.36 0.10 1.00 0.42 0.26 -0.23 0.17
+-0.83 -0.22 0.42 1.00 0.82 0.41 0.02
+-0.86 -0.33 0.26 0.82 1.00 0.52 0.03
+-0.63 -0.75 -0.23 0.41 0.52 1.00 -0.05
+-0.24 -0.30 0.17 0.02 0.03 -0.05 1.00
diff --git a/data/MINERvA/CCQE/proton_Q2QE_nu_data.txt b/data/MINERvA/CC0pi/proton_2014/proton_Q2QE_nu_data.txt
old mode 100755
new mode 100644
similarity index 100%
rename from data/MINERvA/CCQE/proton_Q2QE_nu_data.txt
rename to data/MINERvA/CC0pi/proton_2014/proton_Q2QE_nu_data.txt
diff --git a/data/MINERvA/CC0pi/proton_2014/readme.txt b/data/MINERvA/CC0pi/proton_2014/readme.txt
new file mode 100644
index 0000000..49298bf
--- /dev/null
+++ b/data/MINERvA/CC0pi/proton_2014/readme.txt
@@ -0,0 +1,3 @@
+arxiv.org/abs/1409.4497
+
+muon measurement missing covariance in publication so not included
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_nu_proton.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_nu_proton.cxx
index 312fc99..87fc775 100644
--- a/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_nu_proton.cxx
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_nu_proton.cxx
@@ -1,99 +1,106 @@
// 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
#include
#include "MINERvA_SignalDef.h"
#include "MINERvA_CC0pi_XSec_1DQ2_nu_proton.h"
+// https://arxiv.org/abs/1409.4497
+// DOI: 10.1103/PhysRevD.91.071301
+
//********************************************************************
MINERvA_CC0pi_XSec_1DQ2_nu_proton::MINERvA_CC0pi_XSec_1DQ2_nu_proton(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CC0pi_XSec_1DQ2_nu_proton sample. \n" \
"Target: CH \n" \
- "Flux: MINERvA Forward Horn Current nue + nuebar \n" \
- "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
+ "Flux: MINERvA Forward Horn Current numu \n" \
+ "Signal: Any event with 1 muon, at least one proton in range, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{QE} (GeV^{2})");
fSettings.SetYTitle("d#sigma/dQ^{2} (cm^{2}/GeV^{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MINERvA_CC0pi_XSec_1DQ2_nu_proton");
- fSettings.SetDataInput( FitPar::GetDataBase() + "/MINERvA/CCQE/proton_Q2QE_nu_data.txt" );
- fSettings.SetCovarInput( FitPar::GetDataBase() + "/MINERvA/CCQE/proton_Q2QE_nu_covar.txt" );
+ fSettings.SetDataInput( FitPar::GetDataBase() + "/MINERvA/CC0pi/proton_2014/proton_Q2QE_nu_data.txt" );
+ // Have a shape covariance
+ fIsShape = fSettings.Found("type","SHAPE");
+ std::string covid = FitPar::GetDataBase() + "/MINERvA/CC0pi/proton_2014/proton_Q2QE_nu_covar";
+ if (fIsShape) covid += "_shape";
+ covid += ".txt";
+ fSettings.SetCovarInput(covid);
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38/(fNEvents+0.))/TotalIntegratedFlux();
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCorrelationFromTextFile(fSettings.GetCovarInput() );
// Final setup ---------------------------------------------------
FinaliseMeasurement();
-
};
void MINERvA_CC0pi_XSec_1DQ2_nu_proton::FillEventVariables(FitEvent *event){
// Has NuMuCC1p
if (event->HasISNuMuon() &&
event->HasFSMuon() &&
- event->HasFSProton()){
+ event->HasFSProton()) {
TLorentzVector pnu = event->GetHMISNuMuon()->fP;
TLorentzVector pprot = event->GetHMFSProton()->fP;
TLorentzVector pmu = event->GetHMFSMuon()->fP;
// Q2QE rec from leading proton assuming 34 MeV Eb
double protmax = pprot.E();
double q2qe = FitUtils::ProtonQ2QErec(protmax, 34.);
// Coplanar is angle between muon and proton plane
TVector3 plnprotnu = pprot.Vect().Cross(pnu.Vect());
TVector3 plnmunu = pmu.Vect().Cross(pnu.Vect());
double copl = plnprotnu.Angle(plnmunu);
// Fill X Variables
fXVar = q2qe;
// Save Coplanar into spare y variable
fYVar = copl;
}
return;
};
bool MINERvA_CC0pi_XSec_1DQ2_nu_proton::isSignal(FitEvent *event){
return SignalDef::isCC0pi1p_MINERvA(event, EnuMin*1.E3, EnuMax*1.E3);
};
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_antinu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_antinu.cxx
deleted file mode 100755
index 36d7e17..0000000
--- a/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_antinu.cxx
+++ /dev/null
@@ -1,112 +0,0 @@
-//Adrian Orea
-//I used the file MINERvA_CCinc_XSec_2DEavq3_nu.cxx as a template
-//Also, I am fully aware of the naming typo (should be ptpz), but Everything is already named the same way so...
-
-//Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
-
-/*******************************************************************************
-* This file is part of NUISANCE.
-*
-* NUISANCE is free software: you can redistribute it and/or modify
-* it under the terms of the GNU General Public License as published by
-* the Free Software Foundation, either version 3 of the License, or
-* (at your option) any later version.
-*
-* NUISANCE is distributed in the hope that it will be useful,
-* but WITHOUT ANY WARRANTY; without even the implied warranty of
-* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-* GNU General Public License for more details.
-*
-* You should have received a copy of the GNU General Public License
-* along with NUISANCE. If not, see .
-*******************************************************************************/
-/*
- Author : Adrian Orea
-*/
-
-#include "MINERvA_SignalDef.h"
-#include "MINERvA_CC0pi_XSec_2Dptpx_antinu.h"
-
-//********************************************************************
-MINERvA_CC0pi_XSec_2Dptpx_antinu::MINERvA_CC0pi_XSec_2Dptpx_antinu(nuiskey samplekey) {
-//********************************************************************
-
- // Sample overview ---------------------------------------------------
- std::string descrip = "MINERvA_CC0pi_XSec_2Dptpx_antinu sample. \n" \
- "Target: CH \n" \
- "Flux: MINERvA Medium Energy FHC numu \n" \
- "Signal: CC-0pi \n";
-
- // Setup common settings
- fSettings = LoadSampleSettings(samplekey);
- fSettings.SetDescription(descrip);
- fSettings.SetXTitle("p_{t} (GeV)");
- fSettings.SetYTitle("p_{z} (GeV)");
- fSettings.SetZTitle("d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)");
- fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL");
- fSettings.SetEnuRange(0.0, 100.0);
- fSettings.DefineAllowedTargets("C,H");
-
- // CCQELike plot information
- fSettings.SetTitle("MINERvA_CC0pi_XSec_2Dptpx_antinu");
-
- fSettings.SetDataInput( FitPar::GetDataBase() + "MINERvA/CC0pi/data_2D.txt" );
- fSettings.SetErrorInput( FitPar::GetDataBase() + "MINERvA/CC0pi/error_2D.txt" );
- fSettings.SetCovarInput( FitPar::GetDataBase() + "MINERvA/CC0pi/covar_2D.txt" );
- fSettings.SetMapInput( FitPar::GetDataBase() + "MINERvA/CC0pi/map_2D.txt" );
- fSettings.DefineAllowedSpecies("antinumu");
- FinaliseSampleSettings();
-
- // Scaling Setup ---------------------------------------------------
- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
- fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux();
-
- // Plot Setup -------------------------------------------------------
- Double_t P_t[7] = {0,0.15,0.25,0.4,0.7,1.0,1.5};
- Double_t P_z[12] = {1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,6.0,8.0,10.0,15.0};
-
- CreateDataHistogram(7, P_t, 12, P_z);
- SetDataValuesFromTextFile( fSettings.GetDataInput() );
- ScaleData(1E-41);
-
- SetMapValuesFromText( fSettings.GetMapInput() );
-
- SetDataErrorsFromTextFile( fSettings.GetErrorInput() );
- ScaleDataErrors(1E-41);
-
- //SetCholDecompFromTextFile( fSettings.GetCovarInput() );
- //ScaleCovar(1E-16);
-
- //StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, fMapHist, 1E-38);
-
- // Final setup ---------------------------------------------------
- FinaliseMeasurement();
-
-};
-
-//********************************************************************
-void MINERvA_CC0pi_XSec_2Dptpx_antinu::FillEventVariables(FitEvent *event) {
-//********************************************************************
-
-// Checking to see if there is a Muon
- if (event->NumFSParticle(-13) == 0)
- return;
-
- TLorentzVector Pmu = event->GetHMFSParticle(-13)->fP; //I Added this part
- Double_t px = Pmu.X()/1000;
- Double_t py = Pmu.Y()/1000;
- Double_t pz = Pmu.Z()/1000;
- Double_t pt = sqrt(px*px+py*py);
-
-// Set Hist Variables
- fYVar = pz;
- fXVar = pt;
-
- return;
-};
-
-//********************************************************************
-bool MINERvA_CC0pi_XSec_2Dptpx_antinu::isSignal(FitEvent *event) {
-//********************************************************************
- return (SignalDef::isCC0pi_anti_MINERvAPTPZ(event, -14, EnuMin, EnuMax)/* && SignalDef::HasProtonMomBelowThreshold(event, 120)*/);
-};
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_antinu.h b/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_antinu.h
deleted file mode 100755
index 1cbaac6..0000000
--- a/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_antinu.h
+++ /dev/null
@@ -1,69 +0,0 @@
-//Adrian Orea
-//NEED TO MODIFY
-//I used the file MINERvA_CCinc_XSec_2DEavq3_nu.h as a template
-
-// 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 MINERVA_CC0PI_XSEC_2DPTPX_ANTINU_H_SEEN
-#define MINERVA_CC0PI_XSEC_2DPTPX_ANTINU_H_SEEN
-
-#include "Measurement2D.h"
-
-//********************************************************************
-class MINERvA_CC0pi_XSec_2Dptpx_antinu : public Measurement2D {
-//********************************************************************
-
- public:
-
- // Constructor
- MINERvA_CC0pi_XSec_2Dptpx_antinu(nuiskey samplekey);
-
- // Destructor
- virtual ~MINERvA_CC0pi_XSec_2Dptpx_antinu() {
-
- // Remove all the content histograms *
- // for (int i = 0; i < 9; i++)
- // delete this->fMCHist_content[i];
-
- // delete everything
- /* delete difHist; */
- /* delete evtsignalHist; */
- /* delete fluxsignalHist; */
- /* delete fMapHist; */
- /* delete status; */
- /* delete PDGHistogram; */
-
- /* // Delete MODE histograms */
- /* for (int i = 0; i < 60; i++) */
- /* delete fMCHist_PDG[i]; */
-
- return;
- };
-
- // Required functions
- bool isSignal(FitEvent *nvect);
- void FillEventVariables(FitEvent *event);
-
- protected:
-
- // Cuts
-};
-
-#endif
diff --git a/src/MINERvA/MINERvA_SignalDef.cxx b/src/MINERvA/MINERvA_SignalDef.cxx
index f4f0005..5824f29 100644
--- a/src/MINERvA/MINERvA_SignalDef.cxx
+++ b/src/MINERvA/MINERvA_SignalDef.cxx
@@ -1,468 +1,476 @@
// 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 "SignalDef.h"
#include "FitUtils.h"
#include "MINERvA_SignalDef.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) {
// *********************************
// 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);
// 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;
return true;
};
// Updated MINERvA 2017 Signal using Wexp and no restriction on angle
bool isCC1pip_MINERvA_2017(FitEvent *event, double EnuMin, double EnuMax){
// 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;
// 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) {
// *********************************
// 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.;
// 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
}
if (Wrec > 1800. || Wrec < 0.0) return false;
return true;
};
//********************************************************************
bool isCCQEnumu_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace) {
//********************************************************************
if (!isCCQELike(event, 14, EnuMin, EnuMax)) return false;
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), 34., true);
// If Restricted phase space
if (!fullphasespace && ThetaMu > 0.34906585) return false;
// restrict energy range
if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false;
return true;
};
//********************************************************************
bool isCCQEnumubar_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace) {
//********************************************************************
if (!isCCQELike(event, -14, EnuMin, EnuMax)) return false;
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);
// If Restricted phase space
if (!fullphasespace && ThetaMu > 0.34906585) return false;
// restrict energy range
if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false;
return true;
}
//********************************************************************
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;
// Cut on muon energy < 1.5 GeV
if (pmu.E()/1000.0 < 1.5) return false;
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) {
- // Require numu CC0pi event with a proton above threshold
- bool signal = (isCC0pi(event, 14, enumin, enumax) &&
- HasProtonKEAboveThreshold(event, 110.0));
+ 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;
}
} 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){
//********************************************************************
// 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 ) {
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++;
}
}
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 || abs(pdg) == 111 || abs(pdg) == 130 || abs(pdg) == 310 || abs(pdg) == 311 || abs(pdg) == 313 ) {
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++;
}
}
// 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;
// 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++;
}
}
// Proton momentum cuts
//if (pp.Vect().Mag() < 450 || pp.Vect().Mag() > 1200) return false;
// Proton angle cuts
//if (pp.Vect().Angle(pnu.Vect()) > (M_PI/180.0)*70) return false;
if (nProtonsAboveThreshold == 0) return false;
return true;
};
}