3 #include "CAFAna/Core/Var.h" 5 #include "CAFAna/Core/Binning.h" 42 std::string sFD_ND = ( IsFarDet ==
true ?
"FD" :
"ND" );
48 nom =
"prod_caf_R17-11-14-prod4reco.neutron-respin.c_fd_genie_nonswap_rhc_nova_v08_full_v1";
51 nom =
"prod_caf_R17-11-14-prod4reco.neutron-respin.c_nd_genie_nonswap_rhc_nova_v08_full_v1";
62 const std::vector<std::string>
GENIEStr = {
"All",
"MEC",
"DIS",
"RES",
"COH",
"QE"};
63 const unsigned int nGENIE = GENIECuts.size();
66 std::vector<std::string> cutNames; std::vector<Cut> kDetCuts;
68 cutNames.emplace_back(
"kNoCut"); kDetCuts.emplace_back(
kNoCut);
69 cutNames.emplace_back(
"qual"); kDetCuts.emplace_back(
kNumuQuality);
72 cutNames.emplace_back(
"full"); kDetCuts.emplace_back(
kNumuCutFD2018);
75 cutNames.emplace_back(
"kNoCut"); kDetCuts.emplace_back(
kNoCut);
76 cutNames.emplace_back(
"qual"); kDetCuts.emplace_back(
kNumuQuality);
80 unsigned int nDet = kDetCuts.size();
84 if (!cvmfs_dir) {
std::cerr <<
"Couldn't find UPS dir for numu prediction!" <<
std::endl;
return; }
87 TFile*
inFile = TFile::Open(fdspecfile.c_str());
88 TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny(
"FDSpec2D" );
89 const int NHadEFracQuantiles = 4;
91 HadEFracQuantCuts.push_back(
kNoCut );
92 std::vector<std::string> quantNames = {
"Quant1",
"Quant2",
"Quant3",
"Quant4",
"AllQuant" };
93 const unsigned int nQ = quantNames.size();
98 if (sr->
mc.
nnu == 0)
return -5.;
104 if (sr->
mc.
nnu == 0)
return -5.;
105 if(
kTrueE(sr) == 0)
return -5.;
120 if(sr->
mc.
nnu == 0)
return W;
121 if(sr->
mc.
nu[0].W2<0)
return W;
129 if(sr->
mc.
nnu == 0)
return pos;
130 pos = sr->
mc.
nu[0].vtx.x;
137 if(sr->
mc.
nnu == 0)
return pos;
138 pos = sr->
mc.
nu[0].vtx.y;
145 if(sr->
mc.
nnu == 0)
return pos;
146 pos = sr->
mc.
nu[0].vtx.z;
187 for (
unsigned int quant=0; quant<
nQ; ++quant) {
188 for (
unsigned int gen=0; gen<
nGENIE; ++gen){
189 Cut nom_cut = GENIECuts[gen] && kDetCuts[
det] && HadEFracQuantCuts[quant];
190 std::string nomStr =
" for Cut "+cutNames[
det]+
", "+quantNames[quant] +
", "+GENIEStr[gen];
192 RmT_nom [
det][quant][gen] =
new Spectrum(
"Reco - True (GeV)"+nomStr, RmT_bins, ldr, kRmT, nom_cut,
kNoShift, weight);
193 RmToT_nom [
det][quant][gen] =
new Spectrum(
"(Reco - True)/True"+nomStr, RmToT_bins, ldr, kRmToT, nom_cut,
kNoShift, weight);
194 trueHadE_nom [
det][quant][gen] =
new Spectrum(
"true nu - reco muon"+nomStr, had_bins, ldr,
kTrueE-
kMuE, nom_cut);
195 recoHadE_nom [
det][quant][gen] =
new Spectrum(
"Hadronic Energy"+nomStr, had_bins, ldr,
kHadE, nom_cut);
196 png_nom [
det][quant][gen] =
new Spectrum(
"Number of Prongs"+nomStr, png_bins, ldr, kNprongs, nom_cut);
197 trueNuE_nom [
det][quant][gen] =
new Spectrum(
"kTrueE"+nomStr, mu_bins, ldr,
kTrueE, nom_cut);
198 recoNuE_nom [
det][quant][gen] =
new Spectrum(
"kCCE"+nomStr, mu_bins, ldr,
kCCE, nom_cut);
199 trueMuonE_nom [
det][quant][gen] =
new Spectrum(
"kTrueMuonE"+nomStr, mu_bins, ldr,
kTrueMuonE, nom_cut);
200 recoMuonE_nom [
det][quant][gen] =
new Spectrum(
"kMuE"+nomStr, mu_bins, ldr,
kMuE, nom_cut);
201 remid_nom [
det][quant][gen] =
new Spectrum(
"remid"+nomStr, id_bins, ldr,
SIMPLEVAR(sel.remid.pid), nom_cut);
202 cvnprod3_nom [
det][quant][gen] =
new Spectrum(
"cvnprod3"+nomStr, id_bins, ldr,
SIMPLEVAR(sel.cvnProd3Train.numuid), nom_cut);
203 cvn2017_nom [
det][quant][gen] =
new Spectrum(
"cvn2017"+nomStr, id_bins, ldr,
SIMPLEVAR(sel.cvn2017.numuid), nom_cut);
204 Vx_nom [
det][quant][gen] =
new Spectrum(
"true vertex x"+nomStr, xybins, ldr, kVx, nom_cut);
205 Vy_nom [
det][quant][gen] =
new Spectrum(
"true vertex y"+nomStr, xybins, ldr, kVy, nom_cut);
206 Vz_nom [
det][quant][gen] =
new Spectrum(
"true vertex z"+nomStr, zbins, ldr, kVz, nom_cut);
219 TFile *
OutFile =
new TFile(sOutFile.c_str(),
"RECREATE");
221 for (
unsigned int quant=0; quant<
nQ; ++quant) {
222 for (
unsigned int gen=0; gen<
nGENIE; ++gen) {
223 std::string nomStr = cutNames[
det]+
"_"+quantNames[quant]+
"_"+GENIEStr[gen];
224 RmT_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"RmT_nom_")+TString(nomStr));
225 RmToT_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"RmToT_nom_")+TString(nomStr));
226 trueHadE_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"trueHadE_nom_")+TString(nomStr));
227 recoHadE_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"recoHadE_nom_")+TString(nomStr));
228 png_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"png_nom_")+TString(nomStr));
229 trueNuE_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"trueNuE_nom_")+TString(nomStr));
230 recoNuE_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"recoNuE_nom_")+TString(nomStr));
231 trueMuonE_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"trueMuonE_nom_")+TString(nomStr));
232 recoMuonE_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"recoMuonE_nom_")+TString(nomStr));
233 remid_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"remid_nom_")+TString(nomStr));
234 cvnprod3_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"cvnprod3_nom_")+TString(nomStr));
235 cvn2017_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"cvn2017_nom_")+TString(nomStr));
236 Vx_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"Vx_nom_")+TString(nomStr));
237 Vy_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"Vy_nom_")+TString(nomStr));
238 Vz_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"Vz_nom_")+TString(nomStr));
239 trklen_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"trklen_nom_")+TString(nomStr));
240 scatLL_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"scatLL_nom_")+TString(nomStr));
241 dedxLL_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"dedxLL_nom_")+TString(nomStr));
242 MF_nom [
det][quant][gen] ->
SaveTo(OutFile, TString(
"MF_nom_")+TString(nomStr));
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)
void hyperon_nom_macro(bool IsFarDet=true)
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
void SetSpillCut(const SpillCut &cut)
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.
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
const unsigned int nGENIE
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'...
const SystShifts kNoShift
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;})
const std::vector< std::string > GENIEStr
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...
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
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.