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);
292 const Binning kFDXYBins = Binning::Simple(55,-8.0,8.0);
293 const Binning kFDZBins = Binning::Simple(55,0,60.00);
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(
340 HistAxis(
"Slice N_{Hit}", Binning::Simple(50, 0, 500),
kNHit),
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(
372 HistAxis(
"On-track Calorimetric Hadronic Energy [GeV]", Binning::Simple(50, 0, 0.5),
kNumuHadTrkE),
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(
392 HistAxis(
"Number of Tracks in Slice", Binning::Simple(14, 1, 15),
kNKalman),
393 "reco-nkal",
"Number of tracks in slice. " );
395 variables.emplace_back(
396 HistAxis(
"Number of Hits in Slice", Binning::Simple(50, 0, 500),
kNHit),
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(
460 HistAxis(
"Number of Hits on Primary Track", Binning::Simple(50,0,500),
kTrkNhits),
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(
576 HistAxis(
"Kalman Track Cos #theta_{X}", Binning::Simple(100, -1, 1),
kDirX),
577 "reco-dirX",
"X-direction of muon track. " );
579 variables.emplace_back(
580 HistAxis(
"Kalman Track Cos #theta_{Y}", Binning::Simple(100, -1, 1),
kDirY),
581 "reco-dirY",
"Y-direction of muon track. " );
583 variables.emplace_back(
584 HistAxis(
"Kalman Track Cos #theta_{Z}", Binning::Simple(50, 0, 1),
kDirZ),
585 "reco-dirZ",
"Z-direction of muon track. " );
587 variables.emplace_back(
588 HistAxis(
"Kalman Track Cos #theta_{NuMI}", Binning::Simple(50, 0, 1),
kCosNumi),
589 "reco-cosNumi",
"Beam direction of muon track. " );
591 variables.emplace_back(
592 HistAxis(
"Number of Missing Planes in Primary Track", Binning::Simple(50, 0, 50),
594 "reco-trkNPlaneGap",
"Track N Plane Gap. " );
596 variables.emplace_back(
597 HistAxis(
"Visible Slice Energy Per Hit [GeV]", Binning::Simple(40, 0, 0.04),
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(
609 HistAxis(
"Non QE Neutrino Energy [GeV]", Binning::Simple(50, 0, 5),
kNonQeE),
610 "reco-nonQeE",
"Non QE Energy. " );
612 variablesForQePId.emplace_back(
614 "reco-nonQeHadE",
"Hadronic Non QE Energy. " );
616 variablesForQePId.emplace_back(
617 HistAxis(
"QePId Input: Off-track Energy Ratio", Binning::Simple(100, 0, 1),
kQepidOffE),
618 "reco-QepidOffE",
"Off-track Energy Ratio, QePId. " );
620 variablesForQePId.emplace_back(
621 HistAxis(
"QePId Input: Fractional Energy Difference", Binning::Simple(100, -1, 1),
623 "reco-QepidEDiff",
"Fractional Energy Difference, QePId. " );
625 variablesForQePId.emplace_back(
626 HistAxis(
"QePId Input: Fractional Energy Difference Z-test", Binning::Simple(100, -5, 5),
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(
635 HistAxis(
"QE Neutrino Energy [GeV]", Binning::Simple(50, 0, 5),
kQeE),
636 "reco-qeE",
"QE energy. " );
638 variablesForQePId.emplace_back(
639 HistAxis(
"Hadronic QE Energy [GeV]", Binning::Simple(50, 0, 5),
kQeHadE),
640 "reco-qeHadE",
"Hadronic energy for QE things. " );
642 variablesForQePId.emplace_back(
643 HistAxis(
"QE Angle Neutrino Energy [GeV]", Binning::Simple(50, 0, 5),
kQeAngleE),
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";
const Var kDirY({"trk.nkalman","trk.kalman.dir.fY","sel.remid.bestidx"}, [](const caf::StandardRecord *sr){if(sr->trk.nkalman< 1) return-5.;return sr->trk.kalman[sr->sel.remid.bestidx].dir.Y();})
const Var kSlcExtentY({"slc.boxmax.fY","slc.boxmin.fY"}, [](const caf::StandardRecord *sr){return(sr->slc.boxmax.Y()-sr->slc.boxmin.Y())/100.;})
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;})
const Var kQeHadE({"energy.numusimp.trkqeE","energy.mutrkE","trk.nkalman","sel.remid.bestidx"}, [](const caf::StandardRecord *sr){if(sr->trk.nkalman > 0 &&sr->sel.remid.bestidx!=999) return sr->energy.numusimp.trkqeE-sr->energy.mutrkE[sr->sel.remid.bestidx].E;return-5.f;})
const Var kDirZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 0.f;return sr->trk.kalman.tracks[0].dir.Z();})
_HistAxis< Var > HistAxis
Represent the binning of a Spectrum's x-axis.
const Binning kRemidBinning
Binning for plotting remid attractively.
const Cut kIsDytmanMEC({"mc.nnu","mc.nu.mode"}, [](const caf::StandardRecord *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return(sr->mc.nu[0].mode==caf::kMEC);})
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;})
const Var kTrkNPlaneGap({"trk.kalman.nplanegap","trk.nkalman","sel.remid.bestidx"}, [](const caf::StandardRecord *sr){if(sr->trk.nkalman > 0 &&sr->sel.remid.bestidx!=999) return sr->trk.kalman[sr->sel.remid.bestidx].nplanegap;return(unsigned short) 500;})
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 kSlcMaxX({"slc.boxmax.fX"}, [](const caf::StandardRecord *sr){return sr->slc.boxmax.X()/100.;})
const Var kSlcMaxZ({"slc.boxmax.fZ"}, [](const caf::StandardRecord *sr){return sr->slc.boxmax.Z()/100.;})
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.
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 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 kNonQeHadE({"energy.numusimp.trknonqeE","energy.mutrkE","trk.nkalman","sel.remid.bestidx"}, [](const caf::StandardRecord *sr){if(sr->trk.nkalman > 0 &&sr->sel.remid.bestidx!=999) return sr->energy.numusimp.trknonqeE-sr->energy.mutrkE[sr->sel.remid.bestidx].E;return-5.f;})
const Var kSlcStartTime({"slc.starttime"}, [](const caf::StandardRecord *sr){return sr->slc.starttime/1000.;})
const Var kSlcEndTime({"slc.endtime"}, [](const caf::StandardRecord *sr){return sr->slc.endtime/1000.;})
Track finder for cosmic rays.
const Var kDirX({"trk.nkalman","trk.kalman.dir.fX","sel.remid.bestidx"}, [](const caf::StandardRecord *sr){if(sr->trk.nkalman< 1) return-5.;return sr->trk.kalman[sr->sel.remid.bestidx].dir.X();})
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
const Var kSlcMinY({"slc.boxmin.fY"}, [](const caf::StandardRecord *sr){return sr->slc.boxmin.Y()/100.;})
const Var kSlcMinX({"slc.boxmin.fX"}, [](const caf::StandardRecord *sr){return sr->slc.boxmin.X()/100.;})
const Binning kHadronicEnergyBinning
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;})
const Var kSlcMaxY({"slc.boxmax.fY"}, [](const caf::StandardRecord *sr){return sr->slc.boxmax.Y()/100.;})
const Cut kNumuNCRej([](const caf::SRProxy *sr){return(sr->sel.remid.pid >0.75);})
const Var kSlcMeanTime({"slc.meantime"}, [](const caf::StandardRecord *sr){return sr->slc.meantime/1000.;})
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;})
const Binning kEnergyBinning
const Var kSlcExtentX({"slc.boxmax.fX","slc.boxmin.fX"}, [](const caf::StandardRecord *sr){return(sr->slc.boxmax.X()-sr->slc.boxmin.X())/100.;})
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 kSlcUnCalibNHit({"slc.nhit","slc.ncalhit"}, [](const caf::StandardRecord *sr){return sr->slc.nhit-sr->slc.ncalhit;})
const Var kSlcCalEPerNHit({"slc.nhit","slc.calE"}, [](const caf::StandardRecord *sr){if(sr->slc.nhit > 0) return sr->slc.calE/(1.78 *sr->slc.nhit);return-5.;})
const Var kNumuHadVisE({"energy.numusimp.hadcalE","energy.numusimp.hadtrkE"}, [](const caf::StandardRecord *sr){return sr->energy.numusimp.hadcalE+sr->energy.numusimp.hadtrkE;})
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
const Var kCosNumi({"hdr.det","sel.remid.bestidx","trk.kalman.dir.fX","trk.kalman.dir.fY","trk.kalman.dir.fZ"}, [](const caf::StandardRecord *sr){if(sr->trk.nkalman > 0 &&sr->sel.remid.bestidx!=999){if(sr->hdr.det==1){return sr->trk.kalman[sr->sel.remid.bestidx].dir.Dot(beamDirND);}if(sr->hdr.det==2){return sr->trk.kalman[sr->sel.remid.bestidx].dir.Dot(beamDirFD);}}return-5.;})
const Var kSlcMinZ({"slc.boxmin.fZ"}, [](const caf::StandardRecord *sr){return sr->slc.boxmin.Z()/100.;})
Template for Var and SpillVar.
const Var kSlcExtentZ({"slc.boxmax.fZ","slc.boxmin.fZ"}, [](const caf::StandardRecord *sr){return(sr->slc.boxmax.Z()-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);})