8 #include "CAFAna/Core/HistAxis.h" 26 #include "TDirectory.h" 27 #include "TObjString.h" 52 const double kaonNormalization
54 : fData (loaderData, axis, cut, shiftData, wei),
55 fNCTot (loaderMC, axis, cut &&
kIsNC, shiftMC, wei),
60 fTotal (loaderMC, axis, cut &&
kHasNu, shiftMC, wei),
68 fNumuSelData (loaderData, axisNumuPi, cutNumuPi, shiftData, wei),
69 fNumuSelBkg (loaderMC, axisNumuPi, cutNumuPi && (!
kIsNumuCC ||
kIsAntiNu), shiftMC, wei),
70 fNumuSelCCFromPi (loaderMC, axisNumuPi, cutNumuPi &&
kIsNumuCCFromPiDk, shiftMC, wei),
71 fNumuSelCCFromKa (loaderMC, axisNumuPi, cutNumuPi &&
kIsNumuCCFromKaDk, shiftMC, wei),
75 fNumuUncontainData(
Spectrum::Uninitialized()),
76 fNumuUncontainBkg(
Spectrum::Uninitialized()),
77 fNumuUncontainCCFromPi(
Spectrum::Uninitialized()),
78 fNumuUncontainCCFromKa(
Spectrum::Uninitialized()),
79 fNumuUncontainCCFromOther(
Spectrum::Uninitialized()),
81 fNueCCFromPiEstim(
Spectrum::Uninitialized()),
82 fNueCCFromKaEstim(
Spectrum::Uninitialized()),
86 fIsFixedKaScale(true),
87 fKaonNormalization(kaonNormalization)
98 const double kaonNormalization
102 axis, cut, shiftMC, shiftData, wei,
103 axisNumuPi, cutNumuPi, kaonNormalization)
115 const Cut& cutNumuPi,
119 :
fData (loaderData, axis, cut, shiftData, wei),
120 fNCTot (loaderMC, axis, cut &&
kIsNC, shiftMC, wei),
133 fNumuSelData (loaderData, axisNumuPi, cutNumuPi, shiftData, wei),
163 const Cut& cutNumuPi,
170 axis, cut, shiftMC, shiftData, wei,
171 axisNumuPi, cutNumuPi, axisNumuKa, cutNumuKa)
185 for(
int i = 0;
i < aNue1D.size(); ++
i){
215 int nbinsx= hNue3D->GetNbinsX();
216 int nbinsy= hNue3D->GetNbinsY();
217 int nbinsz= hNue3D->GetNbinsZ();
219 for (
int ee = 1; ee < nbinsz+1; ee++){
220 for (
int pz = 1; pz < nbinsx+1; pz++){
221 for (
int pt = 1;
pt < nbinsy+1;
pt++){
223 double nnu = hNue3D->GetBinContent(pz,
pt,ee);
224 if(nnu>0 && ww>0) hNue3D->SetBinContent(pz,
pt,ee, ww*nnu);
229 std::unique_ptr<TH1D> res(hNue3D->ProjectionZ(
UniqueName().c_str()));
245 double pionScale = 1.;
double kaonScale = 1.;
252 int bin_lo_pi = hpionContained->FindBin(0.75);
253 int bin_hi_pi = hpionContained->FindBin(3);
255 int bin_lo_ka = hkaonUncontain->FindBin(4.5);
256 int bin_hi_ka = hkaonUncontain->FindBin(10);
258 while(epsilon>1
E-6) {
261 kaonCorrect.
Scale (kaonScale);
265 std::unique_ptr<TH1D> hrawContained(rawContained.
ToTH1(pot));
267 pionScale = (hrawContained->Integral(bin_lo_pi,bin_hi_pi))/
268 (hpionContained->Integral(bin_lo_pi,bin_hi_pi));
271 pionCorrect.
Scale(pionScale);
277 std::unique_ptr<TH1D> hrawUncontain(rawUncontain.
ToTH1(pot));
279 kaonScale = (hrawUncontain->Integral(bin_lo_ka,bin_hi_ka))/(hkaonUncontain->Integral(bin_lo_ka,bin_hi_ka));
280 epsilon =
abs(epsilon - kaonScale)/kaonScale;
291 TH2*
tmp = (TH2*)hNumu3D->Project3D(
"yx");
292 auto hWeight2D = (TH2*) tmp->Clone();
296 int nbinsx= hNumu3D->GetNbinsX();
297 int nbinsy= hNumu3D->GetNbinsY();
298 int nbinsz= hNumu3D->GetNbinsZ();
300 for (
int pz = 1; pz < nbinsx+1; pz++){
301 for (
int pt = 1;
pt < nbinsy+1;
pt++){
303 for (
int ee = 1; ee < nbinsz+1; ee++){
304 double nnu = hNumu3D->GetBinContent(pz,
pt,ee);
305 double wei = hratioNumu->GetBinContent(ee);
306 double ene = hratioNumu->GetBinCenter (ee);
307 wsum+= nnu * wei*ene;
311 if(wsum > 0) hWeight2D->SetBinContent(pz,
pt,wsum/sum);
389 TDirectory*
tmp = gDirectory;
391 dir = dir->mkdir(name.c_str());
394 TObjString(
"BENDecomp").Write(
"type");
398 v.Write(
"KaonNormalization");
435 std::unique_ptr<BENDecomp>
438 dir = dir->GetDirectory(name.c_str());
464 if(dir->GetListOfKeys()->Contains(
"KaonNormalization")){
465 auto v = *(
TVectorD*)dir->Get(
"KaonNormalization");
466 ret->fKaonNormalization =
v[0];
467 ret->fIsFixedKaScale=
true;
474 ret->fNumuUncontainCCFromOther = *
Spectrum::LoadFrom(dir,
"numuCC_uncontain_from_other");
477 ret->fNumuPiWeights = 0;
486 <<
". No plots will be saved " <<
std::endl;
491 auto c1 =
new TCanvas(
UniqueName().c_str(),
"c1",800,500);
495 hdata->SetMarkerStyle(kFullCircle);
501 auto hs =
new THStack();
502 for (
auto sp:sp_col){
503 auto hist = sp.first.ToTH1(
pot);
505 hist->SetFillColorAlpha(sp.second,0.3);
506 hs->Add(
hist,
"hist");
511 hdata->GetYaxis()->SetTitle(
"10^{3} Events/ 11.0 #times 10^{20} POT");
512 hdata->GetXaxis()->SetTitle(
"Reconstructed neutrino energy (GeV)");
514 hdata->Draw(
"e1 same");
516 TLegend *
leg =
new TLegend(0.6, 0.55, 0.89,0.89);
517 leg->SetHeader(
"Uncontained");
518 TGraph *
dummy =
new TGraph;
519 dummy->SetLineWidth(3);
520 dummy->SetMarkerStyle(kFullCircle); leg->AddEntry(dummy->Clone(),
"#nu_{#mu} data",
"lpe");
521 dummy->SetMarkerStyle(0);
522 dummy->SetLineColor(kBlack); dummy->SetFillColor(kAzure+4); leg->AddEntry(dummy->Clone(),
"#nu_{#mu} from #pi^{#pm}",
"f");
523 dummy->SetLineColor(kBlack); dummy->SetFillColor(kMagenta-9); leg->AddEntry(dummy->Clone(),
"#nu_{#mu} from K^{#pm}",
"f");
524 dummy->SetLineColor(kBlack); dummy->SetFillColor(
kGreen +2 ); leg->AddEntry(dummy->Clone(),
"#nu_{#mu} from Other",
"f");
525 dummy->SetLineColor(kBlack); dummy->SetFillColor(kGray +1 ); leg->AddEntry(dummy->Clone(),
"Background",
"f");
530 c1->Write(
"numuUncontained");
533 hratioNueKa->GetYaxis()->SetTitle(
"#nu_{e} Ka Estimate / Ka MC");
534 hratioNueKa->GetYaxis()->SetRangeUser(0, 0.2 + hratioNueKa->GetMaximum());
535 hratioNueKa->Draw(
"hist");
536 hratioNueKa->Write(
"nue1Dka");
543 auto c1 =
new TCanvas(
UniqueName().c_str(),
"c1",800,600);
547 hdata2->SetMarkerStyle(kFullCircle);
549 std::vector<std::pair<Spectrum, int> > sp_col2 = {{
fNumuSelBkg,kGray +1},
553 auto hs2 =
new THStack();
554 for (
auto sp:sp_col2){
556 hist2->SetFillColorAlpha(sp.second,0.3);
558 hs2->Add(
hist2,
"hist");
562 hdata2->GetYaxis()->SetTitle(
"10^{3} Events/ 11.0 #times 10^{20} POT");
563 hdata2->GetXaxis()->SetTitle(
"Reconstructed neutrino energy (GeV)");
565 hdata2->Draw(
"e1 same");
569 TLegend *leg2 =
new TLegend(0.6, 0.55, 0.76,0.89);
570 leg2->SetHeader(
"Contained");
571 TGraph *
dummy =
new TGraph;
572 dummy->SetLineWidth(3);
573 dummy->SetMarkerStyle(kFullCircle); leg2->AddEntry(dummy->Clone(),
"#nu_{#mu} data",
"lpe");
574 dummy->SetMarkerStyle(0);
575 dummy->SetLineColor(kBlack); dummy->SetFillColor(kAzure +4 ); leg2->AddEntry(dummy->Clone(),
"#nu_{#mu} from #pi^{#pm}",
"f");
576 dummy->SetLineColor(kBlack); dummy->SetFillColor(kMagenta -9 ); leg2->AddEntry(dummy->Clone(),
"#nu_{#mu} from K^{#pm}",
"f");
577 dummy->SetLineColor(kBlack); dummy->SetFillColor(
kGreen +2); leg2->AddEntry(dummy->Clone(),
"#nu_{#mu} from Other",
"f");
578 dummy->SetLineColor(kBlack); dummy->SetFillColor(kGray +1); leg2->AddEntry(dummy->Clone(),
"Background",
"f");
584 c1->Write(
"numuContained");
588 hratioNumu->GetYaxis()->SetTitle(
"#nu_{#mu} (Data-Other) / (MC-Other)");
590 TLine *
line =
new TLine(0.0,1.0,27.0,1.0);
591 line->SetLineWidth(2);
592 line->SetLineColor(
kRed);
593 hratioNumu->Draw(
"hist");
598 c1->Write(
"numu1Dwgt");
600 fNumuPiWeights->GetXaxis()->SetTitle(
"Ancestor forward momentum p_{z} (GeV/c)");
601 fNumuPiWeights->GetYaxis()->SetTitle(
"Ancestor trans. momentum p_{t} (GeV/c)");
607 c1->Write(
"numu2Dwgt");
610 hnumu3d->GetXaxis()->SetTitle(
"Ancestor forward momentum p_{z} (GeV/c)");
611 hnumu3d->GetYaxis()->SetTitle(
"Ancestor trans. momentum p_{t} (GeV/c)");
612 hnumu3d->GetZaxis()->SetTitle(
"Events");
613 hnumu3d->Project3D(
"yx")->Draw(
"colz");
621 hnue3d->GetXaxis()->SetTitle(
"Ancestor forward momentum p_{z} (GeV/c)");
622 hnue3d->GetYaxis()->SetTitle(
"Ancestor trans. momentum p_{t} (GeV/c)");
623 hnue3d->GetZaxis()->SetTitle(
"Events");
626 hnue3d->Project3D(
"yx")->Draw(
"colz");
627 hnue3d->SetTitle(
"");
635 hratioNuePi->GetYaxis()->SetTitle(
"#nu_{e} Pi Estimate / Pi MC");
636 hratioNuePi->GetYaxis()->SetRangeUser(0.94, 0.2+hratioNuePi->GetMaximum());
637 hratioNuePi->SetTitle(
"");
638 hratioNuePi->Draw(
"hist");
643 c1->Write(
"nue1Dpi");
651 auto c1 =
new TCanvas(
UniqueName().c_str(),
"c1",800,500);
654 auto hs =
new THStack();
655 auto hs2 =
new THStack();
657 std::vector<std::pair<Spectrum, int> > sp_col = {
663 std::vector<std::pair<Spectrum, int> > sp_col2 = {
669 for (
auto sp:sp_col){
670 auto hist = sp.first.ToTH1(thePOT, sp.second);
671 hs->Add(
hist,
"hist");
674 for (
auto sp:sp_col2){
675 auto hist = sp.first.ToTH1(thePOT, sp.second);
676 hist->SetLineStyle(2);
677 hs2->Add(
hist,
"hist");
683 fNue .
ToTH1(thePOT,kPink+9,2)->Draw(
"hist same");
686 gStyle->SetTextSize(0.03);
687 TLegend *
leg =
new TLegend(0.7, 0.6, 0.89,0.89);
688 TGraph *
dummy =
new TGraph;
689 dummy->SetLineWidth(3);
690 dummy->SetLineColor(kPink+9); leg->AddEntry(dummy->Clone(),
"Total Beam #nu_{e}",
"l");
691 dummy->SetLineColor(kAzure+4); leg->AddEntry(dummy->Clone(),
"#nu_{e} from Pi",
"l");
692 dummy->SetLineColor(kMagenta+2); leg->AddEntry(dummy->Clone(),
"#nu_{e} from Ka",
"l");
693 dummy->SetLineColor(kGray+1); leg->AddEntry(dummy->Clone(),
"#nu_{e} from Other",
"l");
694 dummy->SetLineColor(kBlack ); dummy->SetLineStyle(2);
695 leg->AddEntry(dummy->Clone(),
"Uncorrected MC",
"l");
699 c1->Write(
"BENresult");
705 hratioNue->GetYaxis()->SetTitle(
"#nu_{e} Estimate / MC");
706 hratioNue->GetYaxis()->SetRangeUser(0, hratioNue->GetMaximum()+0.2);
707 hratioNue->SetLineColor(kPink+9);
708 hratioNuePi->SetLineColor(kAzure+4);
709 hratioNueKa->SetLineColor(kMagenta+2);
710 hratioNue->Draw(
"hist");
711 hratioNue->SetTitle(
"");
712 hratioNuePi->Draw(
"hist same");
713 hratioNueKa->Draw(
"hist same");
722 c1->Write(
"BENratios");
765 unsigned int tempnhit = sr->
slc.
nhit;
840 std::vector<caf::SRVector3D> temptrackstart;
841 std::vector<caf::SRVector3D> temptrackstop;
870 std::vector<caf::SRVector3D> tempshwlidstart;
871 std::vector<caf::SRVector3D> tempshwlidstop;
872 std::vector<caf::SRVector3D> temppngstart;
873 std::vector<caf::SRCVNParticleResult> tempcvnpart;
874 std::vector<float> tempshwlidcale;
914 sr->
mc.
nu[0].rwgt.ppfx.vuniv.clear();
915 sr->
mc.
nu[0].rwgt.ppfx.nvuniv = 0;
Spectrum fNumuUncontainData
T max(const caf::Proxy< T > &a, T b)
const std::vector< Binning > & GetBinnings() const
Spectrum fNueCCFromPiPtPz
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
void MakeWeightsNumuFromPion() const
const Cut kIsNumuCCNotFromKaDk
Spectrum MC_NCTotalComponent() const override
_HistAxis< Var > HistAxis
Far Detector at Ash River.
A 3-vector with more efficient storage than TVector3.
Cuts and Vars for the 2020 FD DiF Study.
Event ID training variables.
void SavePlotsPi(TDirectory *dir)
Spectrum fNumuUncontainCCFromPi
SRHeader hdr
Header branch: run, subrun, etc.
const Cut kIsNueCCFromPiDk
float ndtrklenact
Near detector – muon track length in active region [cm].
float ndtrkcaltranE
Near detector – muon calorimetric energy in transition plane [GeV].
void ResetCosRejInfo(caf::StandardRecord *sr)
void ClearSecondaryTrackInfo(caf::StandardRecord *sr)
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
unsigned int lastplane
last plane
unsigned int firstplane
first plane
float trknonqeE
Track length non-quasielastic neutrino energy [GeV].
Spectrum MC_AntiNueComponent() const override
Simple record of shifts applied to systematic parameters.
SRVector3D boxmax
Maximum coordinates box containing all the hits [cm].
std::vector< SRFuzzyKProng > png
Vector of 3D prong objects.
SRMichelE me
Michel electron branch.
Proxy for caf::StandardRecord.
Collection of SpectrumLoaders for many configurations.
float ndtrkcalcatE
Near detector – muon calorimetric energy in muon catcher [GeV].
SRTrainingBranch training
Extra training information for prototyping PIDs etc.
void NueEstimateFromPi() const
const Eigen::ArrayXd & GetEigen() const
NB these don't have POT scaling. For expert high performance ops only!
const Cut kIsAntiNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg< 0;})
Is this truly an antineutrino?
const Cut kIsBeamNue(CCFlavSel(12, 12))
Select CC .
Spectrum fNueCCFromKaEstim
unsigned int longestidx
index of longest prong
const HistAxis ktpzAxis("Ancestor target forward momentum tpz (GeV/c)", Binning::Simple(60, 0, 60), kTruetpz)
void CenterTitles(TH1 *histo)
static std::unique_ptr< BENDecomp > LoadFrom(TDirectory *dir, const std::string &name)
Spectrum NCTotalComponent() const override
Spectrum fNumuSelCCFromPi
void Scale(double c)
Multiply this spectrum by a constant c.
Spectrum NCComponent() const override
Representation of a spectrum in any variable, with associated POT.
Spectrum NueComponent() const override
Spectrum fNumuSelCCFromKa
Store BDT variables for the short-baseline oscillation study.
Breakpoint ID (BpfId) output.
static void ReduceForBEN2020Decaf(caf::StandardRecord *sr, const caf::SRProxy *srProxy)
void ResetRVPInfo(caf::StandardRecord *sr)
Spectrum fNumuSelCCFromOther
float trkccE
Track length cc neutrino energy [GeV].
SRKalman kalman
Tracks produced by KalmanTrack.
unsigned int nshwlid
number of shwlid showers - either 0 or number of 3d prongs
SRVector3D boxmin
Minimum coordinates box containing all the hits [cm].
Spectrum MC_NumuComponent() const override
const Cut kIsNumuCCNotFromPiDk
Spectrum AntiNueComponent() const override
float ndhadcalcatE
Near detector – hadronic calorimetric energy NOT on the muon track in muon catcher [GeV]...
float trkqeE
Track length quasielastic neutrino energy [GeV].
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
float hadtrkE
Hadronic calorimetric energy on the muon track[GeV].
const Cut kNue2020NDDecafCut
const Cut kIsNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return!sr->mc.nu[0].iscc;})
Is this a Neutral Current event?
void ClearHoughVertexInfo(caf::StandardRecord *sr)
void ClearMichelTruthInfo(caf::StandardRecord *sr)
float calE
Calorimetric energy of the cluster [GeV].
const Cut kIsNueCCFromKaDk
virtual Spectrum NueEstimate() const
Numu energy estimator output.
void SaveTo(TDirectory *dir, const std::string &name) const override
Spectrum NumuComponent() const override
void SaveTo(TDirectory *dir, const std::string &name) const
SRVector3D vtx
Vertex position in detector coordinates. [cm].
short nnu
Number of neutrinos in nu vector (0 or 1)
float E
Neutrino energy, set to match trkccE [GeV].
An SRSlice contains overarching information for a slice.
Spectrum MC_NueComponent() const override
std::vector< float > Spectrum
unsigned int nhit
number of hits
const Cut kIsNumuCCFromPiDk
void MakeWeightsNumuFromKaon() const
void NueEstimateFromKa() const
SRBpfId bpfid
Output from the BreakPointFitter PID (BPFPIdMaker) package.
Base class for the various types of spectrum loader.
unsigned int ncellsfromedge
minimum number of cells to edge of detector
const Cut kNumuUncontainNDDecafCut
Spectrum fNumuUncontainBkg
The StandardRecord is the primary top-level object in the Common Analysis File trees.
void ResetBPFEnergy(caf::StandardRecord *sr)
float ndhadcaltranE
Near detector – hadronic calorimetric energy NOT on the muon track in transition plane [GeV]...
void SavePlots(TDirectory *dir)
SRNumuEnergy numu
Numu energy estimator.
const HistAxis ktpTAxis("Ancestor target transv. momentum tpT (GeV/c)", Binning::Simple(10, 0, 1), kTruetpT)
Axes used in Beam Nue Decomposition.
void SavePlotsKa(TDirectory *dir)
float ndtrkcalactE
Near detector – muon calorimetric energy in active region [GeV].
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
void ResetLEMInfo(caf::StandardRecord *sr)
const Cut kHasNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return true;})
float ndtrklencat
Near detector – muon track length in muon catcher [cm].
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
Spectrum AntiNumuComponent() const override
Spectrum NCAntiComponent() const override
SRIDBranch sel
Selector (PID) branch.
SRElastic elastic
Single vertex found by Elastic Arms.
Spectrum fNumuSelCCFromPiPtPz
const Cut kIsNueCCNotFromKaDk
assert(nhit_max >=nhit_nbins)
Spectrum fNueCCFromPiEstim
SRSlice slc
Slice branch: nhit, extents, time, etc.
double fKaonNormalization
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
SRFuzzyK fuzzyk
Primary 3D prong object.
float recotrkcchadE
Reconstructed hadronic energy for track cc neutrino energy estimator [GeV].
This module creates Common Analysis Files.
void ResetLIDInfo(caf::StandardRecord *sr)
SRVector3D meanpos
Mean position of hits in slice, weighted by charge [cm].
SRParentBranch parent
True parent branch for matching, e.g. MRCC.
T min(const caf::Proxy< T > &a, T b)
const std::vector< std::string > & GetLabels() const
Prevent histograms being added to the current directory.
const Cut kIsNumuCCFromKaDk
unsigned int ncontplanes
number of continuous planes
SRTrackBranch trk
Track branch: nhit, len, etc.
Spectrum fNumuUncontainCCFromOther
SREnergyBranch energy
Energy estimator branch.
SRTrackBase cosmic
Tracks produced by CosmicTrack.
Spectrum Data_Component() const override
std::string UniqueName()
Return a different string each time, for creating histograms.
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
std::vector< SRTrack > tracks
float hadcalE
Hadronic calorimetric energy NOT on the muon track[GeV].
void ClearMultiNuInfo(caf::StandardRecord *sr)
const Cut kIsNueCCNotFromPiDk
std::vector< SRKalmanTrack > tracks
3D Tracks produced by KalmanTrack
Spectrum MC_AntiNumuComponent() const override
Spectrum fNumuUncontainCCFromKa
const Cut kNumuContainNDDecafCut
SRVertexBranch vtx
Vertex branch: location, time, etc.
SRXnue xnuepid
Output from BDT (XnuePID)