10 #include "CAFAna/Cuts/NumuCuts.h" 11 #include "CAFAna/Core/Binning.h" 13 #include "CAFAna/Core/Var.h" 16 #include "CAFAna/Systs/NumuSysts.h" 17 #include "CAFAna/Vars/NumuVars.h" 31 #include "TMultiGraph.h" 36 #include "TPaveText.h" 40 #include "TAttMarker.h" 57 TLatex*
prelim =
new TLatex(.9, .95,
"NO#nuA Preliminary");
58 prelim->SetTextColor(
kBlue);
60 prelim->SetTextSize(2/30.);
61 prelim->SetTextAlign(32);
84 if(
sr->
mc.
nnu == 0)
return false;
105 const Var kQeHadE({
"energy.numusimp.trkqeE",
"energy.mutrkE",
"trk.nkalman",
106 "sel.remid.bestidx"},
114 const Var kNumuTrackE({
"energy.mutrkE",
"trk.nkalman",
"sel.remid.bestidx"},
122 const Var kNumuHadVisE({
"energy.numusimp.hadcalE",
"energy.numusimp.hadtrkE"},
205 const Var kTrkNPlaneGap({
"trk.kalman.nplanegap",
"trk.nkalman",
"sel.remid.bestidx"},
210 return (
unsigned short) 500;
212 const Var kCosNumi({
"hdr.det",
"sel.remid.bestidx",
"trk.kalman.dir.fX",
"trk.kalman.dir.fY",
213 "trk.kalman.dir.fZ"},
226 const Var kDirX({
"trk.nkalman",
"trk.kalman.dir.fX",
"sel.remid.bestidx"},
229 if(
sr->
trk.nkalman < 1)
return -5.;
233 const Var kDirY({
"trk.nkalman",
"trk.kalman.dir.fY",
"sel.remid.bestidx"},
236 if(
sr->
trk.nkalman < 1)
return -5.;
245 const Var kNonQeHadE({
"energy.numusimp.trknonqeE",
"energy.mutrkE",
"trk.nkalman",
246 "sel.remid.bestidx"},
262 fHist(loader, fAxis.
fObj, fSel.
fObj),
274 std::cout <<
"\nrun : ---- Running CAF NuMU Validation.\n";
279 gStyle->SetMarkerStyle(kFullCircle);
280 TGaxis::SetMaxDigits(3);
284 TH1D* spillHist =
new TH1D(
"reco-spills",
";Detector;Spills", 3, 0, 3);
297 std::vector<Selection> selections;
301 "Selected events pass numu containment, cosmic rejection and slice quality cuts. ");
303 "Selected events pass numu containment, cosmic rejection, slice quality, and ReMId cuts. ");
307 "Selected events pass numu containment and slice quality cuts. ");
309 "Selected events pass numu containment, slice quality, and ReMId cuts. ");
312 std::vector<Selection> selectionsForQePId;
317 "Selected events pass numu containment, cosmic rejection, slice quality, ReMId cuts and have 1 track. ");
321 "Selected events pass numu containment, cosmic rejection, slice quality, ReMId cuts and have 2 tracks. ");
326 "Selected events pass numu containment, slice quality, ReMId cuts and have 1 track. ");
330 "Selected events pass numu containment, slice quality, ReMId cuts and have 2 tracks. ");
333 std::vector<TangibleAxis> variables;
335 variables.emplace_back(
336 HistAxis(
"Reconstructed Neutrino Energy [GeV]", kEnergyBinning,
kCCE),
337 "reco-numuE",
"CC energy estimator. ");
339 variables.emplace_back(
341 "reco-slcNHit",
"Number of hits in slice. " );
343 variables.emplace_back(
344 HistAxis(
"Slice Calorimetric Energy [GeV]", kEnergyBinning,
kCaloE),
345 "reco-calE",
"Calorimetric energy of slice. " );
347 variables.emplace_back(
348 HistAxis(
"Reconstructed Hadronic Energy [GeV]", kHadronicEnergyBinning,
kHadE),
349 "reco-hadE",
"Hadronic enery, i.e. numu energy estimate minus muon track energy. " );
351 variables.emplace_back(
353 "reco-numuTrackE",
"Muon track Energy. " );
355 variables.emplace_back(
357 "reco-hadEPerNHit",
"Average energy per hit in hadronic cluster, i.e. had_E/had_n_hit. " );
359 variables.emplace_back(
361 "reco-trkEPerNHit",
"Average energy per hit on primary track, i.e. trk_E/trk_n_hit. " );
363 variables.emplace_back(
365 "reco-hadNHit",
"Number of hits in hadronic cluster. ");
367 variables.emplace_back(
369 "reco-hadCalE",
"Sum of calibrated energy deposits in hadronic cluster. " );
371 variables.emplace_back(
373 "reco-hadTrkE",
"Sum of calibrated hadronic energy deposits on the muon track. " );
375 variables.emplace_back(
377 "reco-hadVisE",
"Sum of calibrated hadronic energy, both on and off the muon track. " );
380 variables.emplace_back(
382 "reco-maxy",
"Maximum Y position of slice. " );
386 variables.emplace_back(
388 "reco-maxy",
"Maximum Y position of slice. " );
391 variables.emplace_back(
393 "reco-nkal",
"Number of tracks in slice. " );
395 variables.emplace_back(
397 "reco-nhit",
"Number of hits in slice. " );
400 variables.emplace_back(
402 "reco-trkStartX",
"Track start x position. " );
404 variables.emplace_back(
406 "reco-trkStartY",
"Track start y position. " );
408 variables.emplace_back(
410 "reco-trkStartZ",
"Track start z position. " );
412 variables.emplace_back(
414 "reco-trkEndX",
"Track stop x position. " );
416 variables.emplace_back(
418 "reco-trkEndY",
"Track stop y position. " );
420 variables.emplace_back(
422 "reco-trkEndZ",
"Track stop z position. " );
424 variables.emplace_back(
426 "reco-sliceDuration",
"Slice duration. " );
430 variables.emplace_back(
432 "reco-trkStartX",
"Track start x position. " );
434 variables.emplace_back(
436 "reco-trkStartY",
"Track start y position. " );
438 variables.emplace_back(
440 "reco-trkStartZ",
"Track start z position. " );
442 variables.emplace_back(
444 "reco-trkEndX",
"Track stop x position. " );
446 variables.emplace_back(
448 "reco-trkEndY",
"Track stop y position. " );
450 variables.emplace_back(
452 "reco-trkEndZ",
"Track stop z position. " );
454 variables.emplace_back(
456 "reco-sliceDuration",
"Slice duration. " );
459 variables.emplace_back(
461 "reco-trkNhits",
"Number of hits on primary track. ");
463 variables.emplace_back(
465 "reco-trkLength",
"Primary track length. " );
467 variables.emplace_back(
469 "reco-scatLL",
"ReMId scattering log log-likelihood for primary track. " );
471 variables.emplace_back(
473 "reco-dedxLL",
"ReMId dE/dx log log-likelihood for primary track. " );
475 variables.emplace_back(
477 "reco-nonHadPlaneFrac",
"ReMId Non-hadronic plane fraction. " );
479 variables.emplace_back(
481 "reco-remid",
"ReMId kNN score. " );
483 variables.emplace_back(
485 "reco-slcUnCalibNHit",
"Uncalibrated slice nhit. " );
487 variables.emplace_back(
489 "reco-slcMeanTime",
"Slice mean time TNS in microseconds. " );
491 variables.emplace_back(
493 "reco-slcStartTime",
"Slice start time TNS in microseconds. " );
495 variables.emplace_back(
497 "reco-slcEndTime",
"Slice end time TNS in microseconds. " );
500 variables.emplace_back(
502 "reco-slcMinX",
"Slice Min X [m]. " );
504 variables.emplace_back(
506 "reco-slcMinY",
"Slice Min Y [m]. " );
508 variables.emplace_back(
510 "reco-slcMinZ",
"Slice Min Z [m]. " );
512 variables.emplace_back(
514 "reco-slcMaxX",
"Slice Max X [m]. " );
516 variables.emplace_back(
518 "reco-slcMaxY",
"Slice Max Y [m]. " );
520 variables.emplace_back(
522 "reco-slcMaxZ",
"Slice Max Z [m]. " );
524 variables.emplace_back(
526 "reco-slcExtentX",
"Slice Extent X [m]. " );
528 variables.emplace_back(
530 "reco-slcExtentY",
"Slice Extent Y [m]. " );
532 variables.emplace_back(
534 "reco-slcExtentZ",
"Slice Extent Z [m]. " );
538 variables.emplace_back(
540 "reco-slcMinX",
"Slice Min X [m]. " );
542 variables.emplace_back(
544 "reco-slcMinY",
"Slice Min Y [m]. " );
546 variables.emplace_back(
548 "reco-slcMinZ",
"Slice Min Z [m]. " );
550 variables.emplace_back(
552 "reco-slcMaxX",
"Slice Max X [m]. " );
554 variables.emplace_back(
556 "reco-slcMaxY",
"Slice Max Y [m]. " );
558 variables.emplace_back(
560 "reco-slcMaxZ",
"Slice Max Z [m]. " );
562 variables.emplace_back(
564 "reco-slcExtentX",
"Slice Extent X [m]. " );
566 variables.emplace_back(
568 "reco-slcExtentY",
"Slice Extent Y [m]. " );
570 variables.emplace_back(
572 "reco-slcExtentZ",
"Slice Extent Z [m]. " );
575 variables.emplace_back(
577 "reco-dirX",
"X-direction of muon track. " );
579 variables.emplace_back(
581 "reco-dirY",
"Y-direction of muon track. " );
583 variables.emplace_back(
585 "reco-dirZ",
"Z-direction of muon track. " );
587 variables.emplace_back(
589 "reco-cosNumi",
"Beam direction of muon track. " );
591 variables.emplace_back(
594 "reco-trkNPlaneGap",
"Track N Plane Gap. " );
596 variables.emplace_back(
599 "reco-slcCalEPerNHit",
"Slice Energy Per Slice NHit. " );
602 std::vector<TangibleAxis> variablesForQePId;
604 variablesForQePId.emplace_back(
606 "reco-qepid",
"QePId kNN score. " );
608 variablesForQePId.emplace_back(
610 "reco-nonQeE",
"Non QE Energy. " );
612 variablesForQePId.emplace_back(
614 "reco-nonQeHadE",
"Hadronic Non QE Energy. " );
616 variablesForQePId.emplace_back(
618 "reco-QepidOffE",
"Off-track Energy Ratio, QePId. " );
620 variablesForQePId.emplace_back(
623 "reco-QepidEDiff",
"Fractional Energy Difference, QePId. " );
625 variablesForQePId.emplace_back(
628 "reco-QepidEDiffZ",
"Fractional Energy Difference Z-test, QePId. " );
630 variablesForQePId.emplace_back(
632 "reco-QepidDedx",
"dE/dx Ratio, applied only to 2 trk events, QePId. " );
634 variablesForQePId.emplace_back(
636 "reco-qeE",
"QE energy. " );
638 variablesForQePId.emplace_back(
640 "reco-qeHadE",
"Hadronic energy for QE things. " );
642 variablesForQePId.emplace_back(
644 "reco-qeAngleE",
"QE angle formula energy. " );
646 variablesForQePId.emplace_back(
648 "reco-numuTrackE",
"Primary track Energy. " );
650 variablesForQePId.emplace_back(
652 "reco-HadVisE",
"Hadronic Energy . " );
655 std::vector<UsefulHist>
hists;
656 hists.reserve(selections.size() * variables.size() + selectionsForQePId.size()
657 * variablesForQePId.size());
659 for(
const auto& sel:selections){
660 for(
const auto& variable:variables){
661 hists.emplace_back(sel, variable, loader);
665 for(
const auto& sel:selectionsForQePId){
666 for(
const auto& variable:variablesForQePId){
667 hists.emplace_back(sel, variable, loader);
671 std::cout <<
"\nrun : --- run loader.\n";
674 double pot = hists[0].fHist.POT();
677 std::cout <<
"\nrun : --- save output.\n";
679 TFile outputfile(kOutputFileName.c_str(),
"RECREATE");
681 for(
const auto&
hist:hists){
684 temp->Write(tempName.c_str());
689 TH1F *pothist =
new TH1F(
"meta-TotalPOT",
"TotalPOT", 1, 0, 1);
690 TH1F *evthist =
new TH1F(
"meta-TotalEvents",
"TotalEvents", 1, 0, 1);
691 pothist->SetBinContent(1, pot);
692 evthist->SetBinContent(1, spillHist->GetEntries());
693 pothist->Write(
"meta-TotalPOT");
694 evthist->Write(
"meta-TotalEvents");
696 std::cout <<
"\nrun : --- close output files.\n";
Near Detector underground.
const Var kHadNHit([](const caf::SRProxy *sr){unsigned int nought=0;if(sr->trk.kalman.ntracks< 1) return nought;return sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;})
_HistAxis< Var > HistAxis
Far Detector at Ash River.
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
SRHeader hdr
Header branch: run, subrun, etc.
const Var kSlcMeanTime([](const caf::SRProxy *sr){return sr->slc.meantime/1000.;})
const Binning kRemidBinning
Binning for plotting remid attractively.
void caf_numu_validation(std::string kInputFileName, std::string kOutputFileName, bool isFD)
const Var kSlcMaxZ([](const caf::SRProxy *sr){return sr->slc.boxmax.Z()/100.;})
const Var kSlcCalEPerNHit([](const caf::SRProxy *sr){if(sr->slc.nhit > 0) return sr->slc.calE/(1.78 *sr->slc.nhit);return-5.;})
const Var kSlcMaxY([](const caf::SRProxy *sr){return sr->slc.boxmax.Y()/100.;})
SRVector3D boxmax
Maximum coordinates box containing all the hits [cm].
float starttime
start time [ns]
const Var kTrkStartY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Y()/100;})
void Preliminary()
Put NOvA Preliminary on plots.
Tangible< Cut > Selection
const Cut kNumuCosmicRej([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 &&
sr->sel.cosrej.numucontpid2019 > 0.535 && sr->slc.nhit< 400);})
const Var kSliceDuration([](const caf::SRProxy *sr){return(sr->slc.endtime-sr->slc.starttime);})
const Cut kNumuContainFD([](const caf::SRProxy *sr){ std::pair< int, int > planes=calcFirstLastLivePlane(sr->slc.firstplane, std::bitset< 14 >(sr->hdr.dibmask));int planestofront=sr->slc.firstplane-planes.first;int planestoback=planes.second-sr->slc.lastplane;return( sr->slc.ncellsfromedge > 1 &&planestofront > 1 &&planestoback > 1 &&sr->sel.contain.kalfwdcell > 10 &&sr->sel.contain.kalbakcell > 10 &&sr->sel.contain.cosfwdcell > 0 &&sr->sel.contain.cosbakcell > 0);})
const Var kSlcExtentY([](const caf::SRProxy *sr){return(sr->slc.boxmax.Y()-sr->slc.boxmin.Y())/100.;})
void SetSpillCut(const SpillCut &cut)
const Var kDirZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].dir.Z();})
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
const Var kSlcMinX([](const caf::SRProxy *sr){return sr->slc.boxmin.X()/100.;})
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
const Var kHadEPerNHit([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 0.0f;int nHit=sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;if(nHit<=0) return 0.0f;float hadE=sr->energy.numu.hadcalE;return hadE/nHit;})
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
const Var kSlcMaxX([](const caf::SRProxy *sr){return sr->slc.boxmax.X()/100.;})
const Cut kNumuContainND([](const caf::SRProxy *sr){return( sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.ncellsfromedge > 1 &&sr->slc.firstplane > 1
&&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1150
&&( sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&( sr->energy.numu.ndhadcalcatE +sr->energy.numu.ndhadcaltranE)< 0.03
&&sr->sel.contain.kalfwdcellnd > 4 &&sr->sel.contain.kalbakcellnd > 8);})
const Var kDirY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].dir.Y();})
Representation of a spectrum in any variable, with associated POT.
const Var kSlcMinY([](const caf::SRProxy *sr){return sr->slc.boxmin.Y()/100.;})
Track finder for cosmic rays.
float endtime
end time [ns]
SRKalman kalman
Tracks produced by KalmanTrack.
SRVector3D boxmin
Minimum coordinates box containing all the hits [cm].
virtual void AddSpillHistogram(TH1 *h, const SpillVar &var, const SpillCut &cut, const SpillVar &wei=kSpillUnweighted)
Uses include counting the total POT or spills in a run.
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Tangible(const T &obj, const std::string &shortName, const std::string &blurb)
const Binning kHadronicEnergyBinning
UsefulHist(Selection sel, TangibleAxis tanAxis, SpectrumLoader &loader, const Cut &bkg=kNoCut)
const Var kTrkStartZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Z()/100;})
SRRemid remid
Output from RecoMuonID (ReMId) package.
Tangible< HistAxis > TangibleAxis
float calE
Calorimetric energy of the cluster [GeV].
const Binning kEnergyBinning
const Var kNonQeHadE([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999) return float(sr->energy.numu.recotrkcchadE);return-5.f;})
const Var kSlcStartTime([](const caf::SRProxy *sr){return sr->slc.starttime/1000.;})
virtual void Go() override
Load all the registered spectra.
const Cut kNumuNCRej([](const caf::SRProxy *sr){return(sr->sel.remid.pid >0.75);})
short nnu
Number of neutrinos in nu vector (0 or 1)
unsigned int nhit
number of hits
const Var kQeHadE([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999) return sr->energy.numu.trkqeE-sr->energy.numu.recomuonE;return-5.f;})
const Var kTrkStartX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.X()/100;})
const Var kTrkEndZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Z()/100;})
unsigned int ncalhit
number of hits with calibration
The StandardRecord is the primary top-level object in the Common Analysis File trees.
const Var kTrkNPlaneGap([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999) return int(sr->trk.kalman.tracks[0].nplanegap);return 500;})
const Var kNumuTrackE({"energy.mutrkE","trk.nkalman","sel.remid.bestidx"}, [](const caf::StandardRecord *sr){if(sr->trk.nkalman > 0 &&sr->sel.remid.bestidx!=999) return sr->energy.mutrkE[sr->sel.remid.bestidx].E;return-5.f;})
const Var kTrkEndY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Y()/100;})
const Var kTrkEndX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.X()/100;})
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kSlcEndTime([](const caf::SRProxy *sr){return sr->slc.endtime/1000.;})
const Var kDirX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].dir.X();})
SRIDBranch sel
Selector (PID) branch.
assert(nhit_max >=nhit_nbins)
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
SRSlice slc
Slice branch: nhit, extents, time, etc.
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
const Var kSlcExtentX([](const caf::SRProxy *sr){return(sr->slc.boxmax.X()-sr->slc.boxmin.X())/100.;})
SRTrackBranch trk
Track branch: nhit, len, etc.
SREnergyBranch energy
Energy estimator branch.
const Var kSlcExtentZ([](const caf::SRProxy *sr){return(sr->slc.boxmax.Z()-sr->slc.boxmin.Z())/100.;})
const Var kSlcUnCalibNHit([](const caf::SRProxy *sr){return sr->slc.nhit-sr->slc.ncalhit;})
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
const Var kCosNumi([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999){if(sr->hdr.det==1){return sr->trk.kalman.tracks[0].dir.Dot(beamDirND);}if(sr->hdr.det==2){return sr->trk.kalman.tracks[0].dir.Dot(beamDirFD);}}return-5.f;})
const Var kSlcMinZ([](const caf::SRProxy *sr){return sr->slc.boxmin.Z()/100.;})
const Var kTrkNhits([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 65535;return int(sr->trk.kalman.tracks[0].nhit);})
float meantime
mean time, weighted by charge [ns]