11 #include "CAFAna/Core/Binning.h" 20 #include "CAFAna/Core/Var.h" 44 std::string sFD_ND = ( IsFarDet ==
true ?
"FD" :
"ND" );
50 fname =
"prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_fd_genie_nonswap_fhc_nova_v08_full_v1_numu2018";
55 fname =
"prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2018";
62 fname =
"prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_fd_genie_nonswap_rhc_nova_v08_full_v1_numu2018";
67 fname =
"prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_nd_genie_nonswap_rhc_nova_v08_full_v1_numu2018";
81 std::string sOutFile =
"Files_"+sFD_ND+
"_"+sFHC_RHC+
".root";
84 const std::vector<Cut> GENIECuts = {
kNoCut };
85 const std::vector<std::string>
GENIEStr = {
"All" };
88 const unsigned int NumGENIE = GENIECuts.size();
91 std::vector<std::string>
CutNames; std::vector<Cut> kDetCuts;
97 unsigned int nDetCuts = kDetCuts.size();
100 std::string fdspecfile =
"/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles/quantiles__" + sFHC_RHC +
"_full__numu2018.root";
103 TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny(
"FDSpec2D" );
105 const int NHadEFracQuantiles = 4;
108 HadEFracQuantCuts.push_back(
kNoCut );
109 std::vector<std::string>
QuantNames = {
"Quant1",
"Quant2",
"Quant3",
"Quant4",
"AllQuant" };
110 const unsigned int nQuants = QuantNames.size();
116 if (sr->
mc.
nnu == 0)
return -5.;
117 return double(sr->
mc.
nu[0].nneutron);
123 if(sr->
mc.
nnu == 0)
return visE;
126 for (
unsigned int PrimL=0; PrimL < sr->
mc.
nu[0].prim.size(); ++PrimL) {
128 if (sr->
mc.
nu[0].prim[PrimL].pdg != 2112)
continue;
129 visE+=sr->
mc.
nu[0].prim[PrimL].daughterVisEinslc;
135 std::vector<double> visE;
136 if(sr->
mc.
nnu == 0)
return visE;
137 for (
unsigned int PrimL=0; PrimL < sr->
mc.
nu[0].prim.size(); ++PrimL){
138 if (sr->
mc.
nu[0].prim[PrimL].pdg != 2112)
continue;
139 visE.push_back(sr->
mc.
nu[0].prim[PrimL].daughterVisEinslc);
163 if (sr->
mc.
nnu == 0)
return false;
170 Spectrum *NuE_nom [NumGENIE][nDetCuts][nQuants];
171 Spectrum *HadE_nom [NumGENIE][nDetCuts][nQuants];
173 Spectrum *NuE_scale_p [NumGENIE][nDetCuts][nQuants];
174 Spectrum *HadE_scale_p [NumGENIE][nDetCuts][nQuants];
176 Spectrum *NuE_scale_m [NumGENIE][nDetCuts][nQuants];
177 Spectrum *HadE_scale_m [NumGENIE][nDetCuts][nQuants];
179 Spectrum *NNeu [NumGENIE][nDetCuts][nQuants];
180 Spectrum *VisE_sum [NumGENIE][nDetCuts][nQuants];
181 Spectrum *VisE [NumGENIE][nDetCuts][nQuants];
184 for (
unsigned int gen=0; gen<NumGENIE; ++gen) {
185 for (
unsigned int det=0;
det<nDetCuts; ++
det) {
186 for (
unsigned int quant=0; quant<nQuants; ++quant) {
188 Cut GENIECut = kRemoveNoTruth && GENIECuts[gen] && kDetCuts[
det] && HadEFracQuantCuts[quant];
189 std::string MyStr =
" for GENIE Int. "+GENIEStr[gen]+
", Cut "+CutNames[
det]+
", "+QuantNames[quant];
197 HadE_scale_p[gen][
det][quant] =
new Spectrum(
"Reco Had E (Gev), Neutron Scale Up" + MyStr,
Binning::Simple(50, 0, 3), loader,
kHadE, GENIECut,
SystShifts(scale, +1), weight);
200 HadE_scale_m[gen][
det][quant] =
new Spectrum(
"Reco Had E (Gev), Neutron Scale Down" + MyStr,
Binning::Simple(50, 0, 3), loader,
kHadE, GENIECut,
SystShifts(scale, -1), weight);
203 VisE_sum [gen][
det][quant] =
new Spectrum(
"Summed visible energy from neutron daughters (GeV)" + MyStr,
Binning::Simple(50, 0, 3), loader, kVisEsum_Ndau, GENIECut,
kNoShift, weight);
204 VisE [gen][
det][quant] =
new Spectrum(
"Visible energy from neutron daughters (GeV)" + MyStr,
Binning::Simple(50, 0, 3), loader, kVisE_Ndau, GENIECut,
kNoShift, weight);
214 TFile *
OutFile =
new TFile(sOutFile.c_str(),
"RECREATE");
217 for (
unsigned int gen=0; gen<NumGENIE; ++gen) {
218 for (
unsigned int det=0;
det<nDetCuts; ++
det) {
219 for (
unsigned int quant=0; quant<nQuants; ++quant) {
220 std::string MyStr = GENIEStr[gen]+
"_"+CutNames[
det]+
"_"+QuantNames[quant];
223 NuE_nom [gen][
det][quant] ->
SaveTo( OutFile, TString(
"NuE_nom_") + TString(MyStr) ) ;
224 HadE_nom [gen][
det][quant] ->
SaveTo( OutFile, TString(
"HadE_nom_") + TString(MyStr) ) ;
225 NuE_scale_p [gen][
det][quant] ->
SaveTo( OutFile, TString(
"NuE_scale_p_") + TString(MyStr) ) ;
226 HadE_scale_p[gen][
det][quant] ->
SaveTo( OutFile, TString(
"HadE_scale_p_") + TString(MyStr) ) ;
227 NuE_scale_m [gen][
det][quant] ->
SaveTo( OutFile, TString(
"NuE_scale_m_") + TString(MyStr) ) ;
228 HadE_scale_m[gen][
det][quant] ->
SaveTo( OutFile, TString(
"HadE_scale_m_") + TString(MyStr) ) ;
230 NNeu[gen][
det][quant] ->
SaveTo( OutFile, TString(
"NNeu_") + TString(MyStr) ) ;
231 VisE_sum[gen][
det][quant] ->
SaveTo( OutFile, TString(
"VisE_sum_") + TString(MyStr) ) ;
232 VisE[gen][
det][quant] ->
SaveTo( OutFile, TString(
"VisE_") + TString(MyStr) ) ;
caf::Proxy< size_t > npng
caf::Proxy< caf::SRFuzzyK > fuzzyk
Cuts and Vars for the 2020 FD DiF Study.
const Var kNPng2D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-5;return(int) sr->vtx.elastic.fuzzyk.npng2d;})
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);})
Simple record of shifts applied to systematic parameters.
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)
std::string pnfs2xrootd(std::string loc, bool unauth)
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
Encapsulate code to systematically shift a caf::SRProxy.
Representation of a spectrum in any variable, with associated POT.
const Var kNPng3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-5;return(int) sr->vtx.elastic.fuzzyk.npng;})
caf::Proxy< caf::SRElastic > elastic
const Cut kNumuCosmicRej2018([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 && sr->sel.cosrej.numucontpid2019 > 0.53 && sr->slc.nhit< 400 && sr->sel.nuecosrej.pngptp< 0.9 );})
const Cut kNumuContainFD2017
const SpillCut kStandardDQCuts([](const caf::SRSpillProxy *spill){if(spill->dcmedgematchfrac==0 &&spill->fracdcm3hits==0 &&spill->nmissingdcmslg==0) return bool(spill->isgoodspill); if(spill->det==caf::kNEARDET && (spill->fracdcm3hits > 0.45|| spill->nmissingdcms > 0)) return false; if(spill->eventincomplete) return false; if(spill->det==caf::kFARDET && spill->nmissingdcmslg > 0) return false; if(spill->det==caf::kFARDET && !spill->ismc && spill->dcmedgematchfrac<=0.2) return false;return true;})
Cut out events with a noisy detector or with parts missing.
virtual void Go() override
Load all the registered spectra.
const std::string QuantNames[NumQuant]
const Var kNNeutrons([](const caf::SRProxy *sr){float NPar=0.;if(DEBUGGING){std::cout<< "\nLooking at event "<< sr->hdr.evt<< ", there are "<< sr->mc.nu.size()<< ", the 0th neutrino is a "<< sr->mc.nu[0].pdg<< ", with Energy "<< sr->mc.nu[0].E<< ". There were a total of "<< sr->mc.nu[0].prim.size()<< " other primaries associated with that neutrino, "<< sr->mc.nu[0].nneutron<< " were neutrons."<< std::endl;if(kIsDytmanMEC(sr)) std::cout<< " The neutrino was a MEC"<< std::endl;else if(kIsDIS(sr)) std::cout<< " The neutrino was a DIS"<< std::endl;else if(kIsRes(sr)) std::cout<< " The neutrino was a Resonance"<< std::endl;else if(kIsCoh(sr)) std::cout<< " The neutrino was a Coherent"<< std::endl;float TotEn=0.0;for(unsigned int PrimL=0;PrimL< sr->mc.nu[0].prim.size();++PrimL){std::cout<< "\t\tLooking at Prim "<< PrimL<< " of "<< sr->mc.nu[0].prim.size()<< ", it was a "<< sr->mc.nu[0].prim[PrimL].pdg<< ", VisE "<< sr->mc.nu[0].prim[PrimL].visE<< ", VisEInslc "<< sr->mc.nu[0].prim[PrimL].visEinslc<< ", E0 "<< sr->mc.nu[0].prim[PrimL].p.E<< ", enteringE "<< sr->mc.nu[0].prim[PrimL].enteringE<< ", totEscE "<< sr->mc.nu[0].prim[PrimL].totEscE<< std::endl;TotEn+=sr->mc.nu[0].prim[PrimL].totEscE;}std::cout<< "\tThe total uncontained energy was "<< TotEn<< std::endl;NPar=sr->mc.nu[0].nneutron;}return NPar;})
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 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;})
void neutKEsyst(bool IsFHC=false, bool IsFarDet=true)
const Binning kNumuCCEOptimisedBinning
Optimised binning for numuCCE from L. Vinton. See docdb 16332. This was close to 'custC' in that talk...
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...
caf::Proxy< size_t > npng2d
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.
std::vector< std::string > CutNames