3 #include "CAFAna/Core/Binning.h" 7 #include "CAFAna/Core/Var.h" 45 std::string sFD_ND = ( IsFarDet ==
true ?
"FD" :
"ND" );
51 hyp =
"jwolcott_hyperon-ccqe-sim_fd_numubar-nonswap_caf";
57 hyp =
"/pnfs/nova/persistent/users/jwolcott/hyperon-ccqe-sim/nd/numubar-hyperon-ccqe_gen.nd.caf.root";
63 std::vector<std::string> cutNames; std::vector<Cut> kDetCuts;
65 cutNames.emplace_back(
"kNoCut"); kDetCuts.emplace_back(
kNoCut);
66 cutNames.emplace_back(
"qual"); kDetCuts.emplace_back(
kNumuQuality);
69 cutNames.emplace_back(
"full"); kDetCuts.emplace_back(
kNumuCutFD2018);
72 cutNames.emplace_back(
"kNoCut"); kDetCuts.emplace_back(
kNoCut);
73 cutNames.emplace_back(
"qual"); kDetCuts.emplace_back(
kNumuQuality);
77 unsigned int nDet = kDetCuts.size();
81 if (!cvmfs_dir) {
std::cerr <<
"Couldn't find UPS dir for numu prediction!" <<
std::endl;
return; }
84 TFile*
inFile = TFile::Open(fdspecfile.c_str());
85 TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny(
"FDSpec2D" );
86 const int NHadEFracQuantiles = 4;
88 HadEFracQuantCuts.push_back(
kNoCut );
89 std::vector<std::string> quantNames = {
"Quant1",
"Quant2",
"Quant3",
"Quant4",
"AllQuant" };
90 const unsigned int nQ = quantNames.size();
94 if(sr->
mc.
nnu==0)
return false;
95 for(
unsigned int pIdx = 0; pIdx < sr->
mc.nu[0].prim.size(); ++pIdx){
96 if(sr->
mc.
nu[0].prim[pIdx].pdg == 3122)
return true;
102 if(sr->
mc.
nnu==0)
return false;
103 for(
unsigned int pIdx = 0; pIdx < sr->
mc.nu[0].prim.size(); ++pIdx){
104 if(sr->
mc.
nu[0].prim[pIdx].pdg == 3222)
return true;
110 if(sr->
mc.
nnu==0)
return false;
111 for(
unsigned int pIdx = 0; pIdx < sr->
mc.nu[0].prim.size(); ++pIdx){
112 if(sr->
mc.
nu[0].prim[pIdx].pdg == 3212)
return true;
118 if(sr->
mc.
nnu==0)
return false;
119 for(
unsigned int pIdx = 0; pIdx < sr->
mc.nu[0].prim.size(); ++pIdx){
120 if(sr->
mc.
nu[0].prim[pIdx].pdg == 3112)
return true;
125 std::vector<std::string> channelNames; std::vector<Cut> kChannelCuts;
126 channelNames = {
"AllChannels",
"Lambda",
"Sigma+",
"Sigma0",
"Sigma-"};
127 kChannelCuts = {
kNoCut, kIsLambda, kIsSigmaP, kIsSigma0, kIsSigmaM};
128 const unsigned int nCh = channelNames.size();
133 if (sr->
mc.
nnu == 0)
return -5.;
139 if (sr->
mc.
nnu == 0)
return -5.;
140 if(
kTrueE(sr) == 0)
return -5.;
155 if(sr->
mc.
nnu == 0)
return W;
156 if(sr->
mc.
nu[0].W2<0)
return W;
164 if(sr->
mc.
nnu == 0)
return pos;
165 pos = sr->
mc.
nu[0].vtx.x;
172 if(sr->
mc.
nnu == 0)
return pos;
173 pos = sr->
mc.
nu[0].vtx.y;
180 if(sr->
mc.
nnu == 0)
return pos;
181 pos = sr->
mc.
nu[0].vtx.z;
222 for (
unsigned int quant=0; quant<
nQ; ++quant) {
224 Cut hyp_cut = kDetCuts[
det] && HadEFracQuantCuts[quant] && kChannelCuts[
chan];
225 std::string hypStr =
" for Cut "+cutNames[
det]+
", "+quantNames[quant] +
", "+channelNames[
chan];
227 RmT_hyp [
det][quant][
chan] =
new Spectrum(
"Reco - True (GeV)"+hypStr, RmT_bins, ldr, kRmT, hyp_cut);
228 RmToT_hyp [
det][quant][
chan] =
new Spectrum(
"(Reco - True)/True"+hypStr, RmToT_bins, ldr, kRmToT, hyp_cut);
230 recoHadE_hyp [
det][quant][
chan] =
new Spectrum(
"Hadronic Energy"+hypStr, had_bins, ldr,
kHadE, hyp_cut);
231 png_hyp [
det][quant][
chan] =
new Spectrum(
"Number of Prongs"+hypStr, png_bins, ldr, kNprongs, hyp_cut);
235 recoMuonE_hyp [
det][quant][
chan] =
new Spectrum(
"kMuE"+hypStr, mu_bins, ldr,
kMuE, hyp_cut);
237 cvnprod3_hyp [
det][quant][
chan] =
new Spectrum(
"cvnprod3"+hypStr, id_bins, ldr,
SIMPLEVAR(sel.cvnProd3Train.numuid), hyp_cut);
238 cvn2017_hyp [
det][quant][
chan] =
new Spectrum(
"cvn2017"+hypStr, id_bins, ldr,
SIMPLEVAR(sel.cvn2017.numuid), hyp_cut);
239 Vx_hyp [
det][quant][
chan] =
new Spectrum(
"true vertex x"+hypStr, xybins, ldr, kVx, hyp_cut);
240 Vy_hyp [
det][quant][
chan] =
new Spectrum(
"true vertex y"+hypStr, xybins, ldr, kVy, hyp_cut);
241 Vz_hyp [
det][quant][
chan] =
new Spectrum(
"true vertex z"+hypStr, zbins, ldr, kVz, hyp_cut);
254 TFile *
OutFile =
new TFile(sOutFile.c_str(),
"RECREATE");
256 for (
unsigned int quant=0; quant<
nQ; ++quant) {
279 RmT_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"RmT_hyp_")+TString(hypStr));
280 RmToT_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"RmToT_hyp_")+TString(hypStr));
281 trueHadE_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"trueHadE_hyp_")+TString(hypStr));
282 recoHadE_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"recoHadE_hyp_")+TString(hypStr));
283 png_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"png_hyp_")+TString(hypStr));
284 trueNuE_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"trueNuE_hyp_")+TString(hypStr));
285 recoNuE_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"recoNuE_hyp_")+TString(hypStr));
286 trueMuonE_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"trueMuonE_hyp_")+TString(hypStr));
287 recoMuonE_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"recoMuonE_hyp_")+TString(hypStr));
288 remid_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"remid_hyp_")+TString(hypStr));
289 cvnprod3_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"cvnprod3_hyp_")+TString(hypStr));
290 cvn2017_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"cvn2017_hyp_")+TString(hypStr));
291 Vx_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"Vx_hyp_")+TString(hypStr));
292 Vy_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"Vy_hyp_")+TString(hypStr));
293 Vz_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"Vz_hyp_")+TString(hypStr));
294 trklen_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"trklen_hyp_")+TString(hypStr));
295 scatLL_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"scatLL_hyp_")+TString(hypStr));
296 dedxLL_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"dedxLL_hyp_")+TString(hypStr));
297 MF_hyp [
det][quant][
chan] ->
SaveTo(OutFile, TString(
"MF_hyp_")+TString(hypStr));
caf::Proxy< size_t > npng
caf::Proxy< caf::SRFuzzyK > fuzzyk
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Cut kNumuContainND2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;}if(sr->trk.kalman.ntracks< 1) return false;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(i==sr->trk.kalman.idxremid) continue;else if(sr->trk.kalman.tracks[i].start.Z() > 1275||sr->trk.kalman.tracks[i].stop.Z() > 1275) return false;}return(sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.firstplane > 1
&&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1100
&&(sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&sr->sel.contain.kalfwdcellnd > 5 &&sr->sel.contain.kalbakcellnd > 10);})
Proxy for caf::StandardRecord.
caf::Proxy< std::vector< caf::SRNeutrino > > nu
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
caf::Proxy< short int > nnu
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Representation of a spectrum in any variable, with associated POT.
const double kAna2018RHCPOT
caf::Proxy< caf::SRElastic > elastic
const Var kTrueE([](const caf::SRProxy *sr){assert(sr->mc.nnu==1);return sr->mc.nu[0].E;})
const Cut kNumuContainFD2017
std::string getenv(std::string const &name)
virtual void Go() override
Load all the registered spectra.
std::vector< float > Spectrum
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to 'custC'...
caf::Proxy< caf::SRTruthBranch > mc
caf::Proxy< bool > IsValid
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kTrueMuonE([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return 0.f;if(sr->mc.nu[0].prim.empty()) return 0.f;if(std::abs(sr->mc.nu[0].prim[0].pdg)!=13) return 0.f;return float(sr->mc.nu[0].prim[0].p.E);})
const Cut kNumuPID2018([](const caf::SRProxy *sr){std::cout<< "ERROR::kNumuPID2018, cutting on both cvnProd3Train and cvn2017."<< " Neither branch exists anymore. Returning False."<< std::endl;abort();return false;})
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
caf::Proxy< caf::SRVertexBranch > vtx
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
void hyperon_macro(bool IsFarDet=true)