29 #include "Utilities/AssociationUtil.h" 30 #include "Utilities/func/MathUtil.h" 34 #include "NovaDAQConventions/DAQConventions.h" 190 mf::LogWarning(
"NumuEAna") <<
"attempting to run on detector that is not " 208 for(
size_t h = 0;
h < hithdl->size(); ++
h){
220 if(slicecol->empty()){
226 for(
unsigned int i = 0;
i<slicecol->size();++
i){
233 mf::LogWarning(
"NumuEAna") <<
"attempting to run MC truth check on " 234 <<
"things without truth, bail";
246 for(
size_t s = 0;
s < slicelist.
size(); ++
s){
251 if (!sliceEnergy)
continue;
254 if (!sliceCosrej)
continue;
269 if(dibmask == ((1<<16)-1)) {
276 unsigned int dibfirst = 0;
277 unsigned int diblast = 0;
280 iD=0;
while (!((dibmask>>iD)&1)) iD++; dibfirst=iD+1;
281 iD=0;
while (dibmask>>iD) iD++; diblast=iD;
284 unsigned int cellsfromedge = 99999999;
291 for(
unsigned int hitIdx = 0; hitIdx < slicelist.
at(
s)->NCell(); ++hitIdx)
294 const int planeNum = chit->
Plane();
345 sliceHits = slicelist[
s]->AllCells();
348 if (funcReturn.empty())
continue;
354 const simb::MCNeutrino& test_neutrino = funcReturn[0].neutrinoInt->GetNeutrino();
364 std::set<int> neutrinoTrackIDs;
365 for (
unsigned int i = 0;
i < particleList.size(); ++
i){
366 neutrinoTrackIDs.insert(particleList[
i]->TrackId());
369 std::set<int>::iterator itr = neutrinoTrackIDs.begin();
378 while( itr != neutrinoTrackIDs.end() ){
384 if ((PDG != 13)&&(PDG != -13)){
391 if (mother->
TrackId() != id_int){
407 double closeToCentre = 50.0;
411 for(
unsigned int q = 0; q < flshits.size(); ++q){
412 if((flshits[q].GetPlaneID() == 191) && (
std::abs(flshits[q].GetPDG())==13)){
413 if(
fabs(flshits[q].GetEntryY()) < closeToCentre){
414 catcherE = flshits[q].GetEntryEnergy() + 0.105658;
415 closeToCentre =
fabs(flshits[q].GetEntryY());
417 if(
fabs(flshits[q].GetExitY()) < closeToCentre){
418 catcherE = flshits[q].GetExitEnergy() + 0.105658;
419 closeToCentre =
fabs(flshits[q].GetExitY());
433 const std::vector< art::Ptr<rb::Track> > sliceTracks = sliceToTracks.at(
s);
434 if(sliceTracks.empty())
continue;
437 const std::vector< art::Ptr<rb::Track> > sliceCosmicTracks = sliceToCosmicTracks.at(
s);
438 nCosmic = sliceCosmicTracks.size();
455 if(trackEnergy.isValid()){
461 std::set<int> muonTrackID;
464 bestTrackHits = sliceTracks[
bestTrack]->AllCells();
516 if(
start_z_reco < 1150 && (end_z_reco < 1275 || sliceCosrej->KalYPosAtTrans() < 55)){
569 fNumuETree = tfs->
make<TTree>(
"NumuETree",
"Numu analysis Energy tree");
640 fCalCCEnergy = tfs->
make<TH1D>(
"CalCCEnergy",
";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
641 fTrkQEEnergy = tfs->
make<TH1D>(
"TrkQEEnergy",
";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
642 fTrkNonQEEnergy = tfs->
make<TH1D>(
"TrkNonQEEnergy",
";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
643 fAngleQEEnergy = tfs->
make<TH1D>(
"AngleQEEnergy",
";Reco Neutrino Energy (GeV);Slices",200, -0.1, 10);
645 fAngleQEError = tfs->
make<TH1D>(
"AngleQEError",
";Error on Neutrino Energy;Slices",200, -0.1, 5.0);
647 fCalCCEnergyDiff = tfs->
make<TH1D>(
"CalCCEnergyDiff",
";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
648 fTrkQEEnergyDiff = tfs->
make<TH1D>(
"TrkQEEnergyDiff",
";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
649 fTrkNonQEEnergyDiff = tfs->
make<TH1D>(
"TrkNonQEEnergyDiff",
";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
650 fAngleQEEnergyDiff = tfs->
make<TH1D>(
"AngleQEEnergyDiff",
";(Reco - Neutrino Energy)/Neutrino Energy;Slices",500, -3.0, 3.0);
652 fCalCCEnergy2D = tfs->
make<TH2D>(
"CalCCEnergy2D",
";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
653 fTrkQEEnergy2D = tfs->
make<TH2D>(
"TrkQEEnergy2D",
";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
654 fTrkNonQEEnergy2D = tfs->
make<TH2D>(
"TrkNonQEEnergy2D",
";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
655 fAngleQEEnergy2D = tfs->
make<TH2D>(
"AngleQEEnergy2D",
";Reco Energy (GeV);Neutrino Energy (GeV)",200, 0., 20, 200, 0., 20);
657 fHadronCalEnergy = tfs->
make<TH1D>(
"HadronCalEnergy",
";Reco Hadron Energy (GeV);Slices",200, -0.1, 10);
658 fHadronTrkEnergy = tfs->
make<TH1D>(
"HadronTrkEnergy",
";Track Hadron Energy (GeV);Slices",200, -0.1, 5);
660 fNDTrkLenActive = tfs->
make<TH1D>(
"NDTrkLenActive",
";Track Length (cm) in Active Region;Slices",200, 0., 1400);
661 fNDTrkLenCatcher = tfs->
make<TH1D>(
"NDTrkLenCatcher",
";Track Length (cm) in Catcher Region;Slices",200, 0., 800);
662 fNDTrkCalActive = tfs->
make<TH1D>(
"NDTrkCalActive",
";Track Energy (GeV) in Active Region;Slices",200, 0., 10);
663 fNDTrkCalTransition = tfs->
make<TH1D>(
"NDTrkCalTransition",
";Track Energy (GeV) in Tran. Plane;Slices",200, 0., 0.5);
664 fNDTrkCalCatcher = tfs->
make<TH1D>(
"NDTrkCalCatcher",
";Track Energy (GeV) in Catcher Region;Slices",200, 0., 2);
665 fNDHadCalActive = tfs->
make<TH1D>(
"NDHadCalActive",
";Hadron Energy (GeV) in Active Region;Slices",200, 0., 10);
666 fNDHadCalTransition = tfs->
make<TH1D>(
"NDHadCalTransition",
";Hadron Energy (GeV) in Tran. Plane;Slices",200, 0., 0.5);
667 fNDHadCalCatcher = tfs->
make<TH1D>(
"NDHadCalCatcher",
";Hadron Energy (GeV) in Catcher Region;Slices",200, 0., 2);
695 mf::LogWarning(
"RecoCheckAna") <<
"Particle has no trajectory points";
703 for(
unsigned int n = 0;
n < N-1; ++
n){
706 if(p0.Z() < 0 || p1.Z() < 0 ||
719 dist += (p1-
p0).
Mag();
double E(const int i=0) const
std::string fPIDModuleLabel
label for module writing pid objects
std::string fCosmicTrackModuleLabel
label for module making cosmic tracks
unsigned int NumberTrajectoryPoints() const
SubRunNumber_t subRun() const
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
TH2D * fCalCCEnergy2D
histogram of calorimetric neutrino E and true neutrino E (for pure pop.)
Energy estimators for CC events.
TH1D * fNDHadCalCatcher
histogram of ND hadron calorimetric energy in the muon catcher region
TTree * fNumuETree
n-tuple used for further analyses (pure population)
fvar< T > fabs(const fvar< T > &x)
TH1D * fHadronCalEnergy
histogram of hadronic calorimetric energy NOT on the muon track
const sim::Particle * TrackIDToMotherParticle(int const &id) const
bool contained
does it pass numu containment cuts?
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
TH1D * fNDTrkLenActive
histogram of ND muon track length in the active region
std::string fEnergyModuleLabel
label for module writing energy objects
const simb::MCParticle & Nu() const
virtual void analyze(art::Event const &e)
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
int qepidntrk
Number of tracks used by QePId.
Vertical planes which measure X.
std::string fTrackModuleLabel
label for module making the tracks
unsigned int Ncells() const
Number of cells in this plane.
art::ServiceHandle< nova::dbi::RunHistoryService > rh
TH1D * fHadronTrkEnergy
histogram of hadronic calorimetric energy on the muon track
double TrueMuonE
True energy of the muon.
TH1D * fNDTrkCalCatcher
histogram of ND muon calorimetric energy in the muon catcher region
TH1D * fCalCCEnergyDiff
histogram of the diff. between calorimetric neutrino E and true neutrino E (for pure pop...
bool trackMatchesMuon
track is muon-like
double cosrejnumucontpid
Cosmic rejection PID for contained events from CosRej.
charged current quasi-elastic
const PlaneGeo * Plane(unsigned int i) const
TH1D * fAngleQEError
histogram of error on the neutrino energy for QE population using angle formula
double RecoMuonE
Reco energy of the muon.
DEFINE_ART_MODULE(TestTMapFile)
TH1D * fNDTrkCalTransition
histogram of ND muon calorimetric energy in the transition region
double bestPID
Highest RemID.
TH1D * fTrkNonQEEnergy
histogram of muon track length neutrino energy for NonQE super pure population
Horizontal planes which measure Y.
std::string fHitModuleLabel
lable for module making the hits
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
double catcherE
Energy before entering the muon catcher (ND)
std::string fSlicerModuleLabel
label for module making the slices
NumuEAna(fhicl::ParameterSet const &p)
unsigned short Cell() const
TH1D * fNDHadCalActive
histogram of ND hadron calorimetric energy in the active region
int InteractionType() const
TH1D * fNDTrkLenCatcher
histogram of ND muon track length in the muon catcher region
Far Detector at Ash River, MN.
bool TrueQEInt
is the interaction QE?
TH2D * fTrkQEEnergy2D
histogram of muon track length neutrino E and true neutrino E (for pure QE pop.)
double cosrejanglekal
Cosine of angle of best ReMId Kalman Track from CosRej.
void push_back(Ptr< U > const &p)
double TrueNuE
True energy of the neutrino.
TH1D * fTrkQEEnergy
histogram of muon track length neutrino energy for QE super pure population
std::string fQepidLabel
label for module writing qe pid objects
double TrueTrackLength
True tracklength of the muon.
virtual void reconfigure(fhicl::ParameterSet const &p)
unsigned int mincellsfromedge
T get(std::string const &key) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Near Detector in the NuMI cavern.
TH1D * fAngleQEEnergy
histogram of neutrino energy for QE population using angle formula
TH1D * fTrkQEEnergyDiff
histogram of the diff. between muon track length neutrino E and true neutrino E (for pure QE pop...
EventNumber_t event() const
reference at(size_type n)
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
double HitCollectionPurity(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, std::map< int, double > *purMap=0, std::map< int, int > *parents=0, bool energyPur=false) const
Returns the fraction of hits in a collection that come from the specified Geant4 track ids...
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
TH2D * fTrkNonQEEnergy2D
histogram of muon track length neutrino E and true neutrino E (for pure NonQE pop.)
TH1D * fNDTrkCalActive
histogram of ND muon calorimetric energy in the active region
T * make(ARGS...args) const
double NuPurity
Neutrino selection purity.
bool isCC
is the interaction CC?
double DetHalfWidth() const
double qepidPID
PID value of output kNN from QePId.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
double pythag(double x, double y)
2D Euclidean distance
double fMaxTrkLen
length of detector from corner to corner
TH1D * fNDHadCalTransition
histogram of ND hadron calorimetric energy in the transition region
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
double trackLength
Reconstructed tracklength.
TH1D * fAngleQEEnergyDiff
histogram of the diff. between neutrino E using angle formula and true neutrino E (for pure QE pop...
bool PassMuonCuts(int trackID, art::PtrVector< rb::CellHit > const &hits) const
Tool for NumuEAna that allows one to see if primary muon (or any track ID you feed the function) cont...
T min(const caf::Proxy< T > &a, T b)
double NuEfficiency
Neutrino selection efficiency.
int GoodDiBlockMask(int subrun=-1, bool reload=false)
std::string fCosRejLabel
label for module writing cosmic rejection
Event generator information.
double TotalLengthInDetector(const sim::Particle *part, geo::View_t view)
virtual void beginRun(art::Run const &r)
TH1D * fTrkNonQEEnergyDiff
histogram of the diff. between muon track length neutrino E and true neutrino E (for pure NonQE pop...
Encapsulate the geometry of one entire detector (near, far, ndos)
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
TH1D * fCalCCEnergy
histogram of calorimetric neutrino energy for general CC population
TH2D * fAngleQEEnergy2D
histogram of neutrino E using angle formula and true neutrino E (for pure QE pop.) ...
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const