diff --git a/CMakeLists.txt b/CMakeLists.txt index 4db7f7000..3bb1590b7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -32,7 +32,7 @@ find_package( ubcv REQUIRED EXPORT ) find_package( larana REQUIRED EXPORT ) find_package( larpandora REQUIRED EXPORT ) find_package( swtrigger REQUIRED EXPORT ) -find_package(Eigen3 REQUIRED) +find_package( Eigen3 REQUIRED ) # macros for dictionary and simple_plugin include(ArtDictionary) diff --git a/ubana/searchingfornues/Selection/AnalysisTools/BlipAnalysis_tool.cc b/ubana/searchingfornues/Selection/AnalysisTools/BlipAnalysis_tool.cc index 685f97ba8..d22c61b0b 100644 --- a/ubana/searchingfornues/Selection/AnalysisTools/BlipAnalysis_tool.cc +++ b/ubana/searchingfornues/Selection/AnalysisTools/BlipAnalysis_tool.cc @@ -1,6 +1,8 @@ #ifndef ANALYSIS_BLIPANALYSIS_CXX #define ANALYSIS_BLIPANALYSIS_CXX +#include "ubreco/BlipReco/RNN/BlipRNNAlg.h" // gives error if added before other includes + #include #include "AnalysisToolBase.h" @@ -14,7 +16,9 @@ #include "ubana/searchingfornues/Selection/CommonDefs/SCECorrections.h" #include "ubana/searchingfornues/Selection/CommonDefs/Geometry.h" #include "lardataobj/RecoBase/SpacePoint.h" -#include "ubobj/Blip/DataTypes.h" + +// blip reco stuff +#include "ubreco/BlipReco/Alg/BlipRecoAlg.h" namespace analysis { @@ -39,6 +43,10 @@ namespace analysis // Original tool created by Afroditi Papadopoulou & Burke Irwin (burke.irwin7@gmail.com) on 01/24/2023. // Updated for MCC10 by Will Foreman (wforeman.phys@gmail.com) on 01/28/2025. // +// Updates (Will Foreman, 02/19/2026) +// - Option to re-run blip reconstruction instead of reading saved blips +// - Inclusion of Liani Silva's proton blip direction RNN +// //////////////////////////////////////////////////////////////////////// class BlipAnalysis : public AnalysisToolBase @@ -62,7 +70,7 @@ class BlipAnalysis : public AnalysisToolBase void analyzeSlice(art::Event const &e, std::vector &slice_pfp_v, bool fData, bool selected) override; /// @brief Method for adding a blip's info to the branch vectors - void addTheBlip(art::Ptr blip); + void addTheBlip(blipobj::Blip const &blip); /// @brief set branches for TTree void setBranches(TTree *_tree) override; @@ -77,6 +85,13 @@ class BlipAnalysis : public AnalysisToolBase float truncate(float input, double base); private: + + + blip::BlipRecoAlg *fBlipAlg; + blip::BlipRNNAlg *fBlipRNN; + TVector3 _blipdir; + bool _blipdir_isValid; + bool _RNNModelLoaded = false; // --- fcl parameters --- art::InputTag fBlipProducer;// blip collection producer @@ -85,9 +100,13 @@ class BlipAnalysis : public AnalysisToolBase bool fSaveSCECorrEnergy; // option to save SCE- and lifetime-corrected energy bool fSaveOnlyNuEvts; // only save blips for neutrino PFP events bool fLiteMode; // save bare minimum of variables + // --- new params 02/2026 --- + bool fRerunBlipReco; // Re-run blip reco and use new output + bool fEnableRNN; // Use blip-direction predictor (TorchLib) // --- event identifiers --- int _run, _sub, _evt; + float _elifetime; // --- 3D blip information --- int _nblips; // number of blips found in this event @@ -96,6 +115,8 @@ class BlipAnalysis : public AnalysisToolBase std::vector _blip_nplanes; // number of planes matched (2 or 3) std::vector _blip_proxtrkdist;// distance to closest track std::vector _blip_proxtrkid; // track ID of closest track + std::vector _blip_trkid; // track ID of track that blips shares hits with + std::vector _blip_trkidfrac; // - fraction of blip's hits that are shared w/track std::vector _blip_touchtrk; // is blip touching a track (ie, delta-like?) std::vector _blip_touchtrkid; // track ID of touched track std::vector _blip_bydeadwire; // is blip adjacent to a dead channel on coll plane? @@ -112,6 +133,11 @@ class BlipAnalysis : public AnalysisToolBase std::vector _blip_true_g4id; // truth-matched G4 track ID std::vector _blip_true_pdg; // truth-matched PDG std::vector _blip_true_energy;// true energy deposited + std::vector _blip_primary_ancestor_pdg; // primary blip ancestor PDG + std::vector _blip_primary_ancestor_g4id; // primary ancestor G4 TrackID + std::vector _blip_true_dir_x; // Initial momentum direction of truth-matched particle + std::vector _blip_true_dir_y; // Initial momentum direction of truth-matched particle + std::vector _blip_true_dir_z; // Initial momentum direction of truth-matched particle // plane-specific information std::vector _blip_pl0_nwires; // number of wires on this plane std::vector _blip_pl1_nwires; // number of wires on this plane @@ -125,6 +151,22 @@ class BlipAnalysis : public AnalysisToolBase std::vector _blip_pl0_centerwire; // central wire number std::vector _blip_pl1_centerwire; // central wire number std::vector _blip_pl2_centerwire; // central wire number + std::vector _blip_rnn_dir_isValid; // Is RNN output valid? + std::vector _blip_rnn_dir_x; // RNN-predicted blip direction (protons) + std::vector _blip_rnn_dir_y; // RNN-predicted blip direction (protons) + std::vector _blip_rnn_dir_z; // RNN-predicted blip direction (protons) + + std::vector _blip_true_category; // Help categorize origin of blip + // -9 = no truth match (data/overlay) + // 0 = truth-matched, but not falling in category + // 1 = primary (n,1p) + // 2 = primary (n,Np) + // 3 = secondary (n,1p) + // 4 = secondary (n,Np) + // 5 = primary (n,gamma) + // 6 = secondary (n,gamma) + // 7 = ncapture gamma + // 8 = mu capture gamma }; @@ -132,12 +174,24 @@ class BlipAnalysis : public AnalysisToolBase //---------------------------------------------------------------------------- BlipAnalysis::BlipAnalysis(const fhicl::ParameterSet &p) { + fRerunBlipReco = p.get ("RerunBlipReco", false); fBlipProducer = p.get("BlipProducer", "blipreco"); fSaveOnlyNuEvts = p.get ("SaveOnlyNuEvts", true); fNuBlipRadius = p.get ("NuBlipRadius", 300.); fSaveSCECorrLoc = p.get ("SaveSCECorrLocation", true); fSaveSCECorrEnergy = p.get ("SaveSCECorrEnergy", true); fLiteMode = p.get ("LiteMode", false); + fEnableRNN = p.get ("EnableRNN", false); + + _blipdir_isValid = false; + if( fRerunBlipReco ) { + fBlipAlg = new blip::BlipRecoAlg( p.get("BlipAlg") ); + if( fEnableRNN ){ + fBlipRNN = new blip::BlipRNNAlg( p.get("BlipRNN") ); + _RNNModelLoaded = fBlipRNN->model_loaded; + } + } + } @@ -153,24 +207,22 @@ float BlipAnalysis::truncate(float input, double base){ } //---------------------------------------------------------------------------- -void BlipAnalysis::addTheBlip(art::Ptr blip){ - +void BlipAnalysis::addTheBlip( blipobj::Blip const &blip) { // get reconstructed position - //TVector3 loc; - TVector3 loc = (fSaveSCECorrLoc) ? blip->PositionSCE : blip->Position; + TVector3 loc = (fSaveSCECorrLoc) ? blip.PositionSCE : blip.Position; // get reconstructed charge and energy float energy = -9; float charge = -9; if( fSaveSCECorrEnergy ) { - energy = blip->Energy; - charge = blip->Charge; + energy = blip.Energy; + charge = blip.Charge; } else { - energy = blip->EnergyCorr; - charge = blip->ChargeCorr; + energy = blip.EnergyCorr; + charge = blip.ChargeCorr; } - float energyTrue = blip->truth.Energy; - float size = sqrt( pow(blip->dX,2) + pow(blip->dYZ,2) ); + float energyTrue = blip.truth.Energy; + float size = sqrt( pow(blip.dX,2) + pow(blip.dYZ,2) ); float x = loc.X(); float y = loc.Y(); float z = loc.Z(); @@ -179,10 +231,10 @@ void BlipAnalysis::addTheBlip(art::Ptr blip){ int nwirestot = 0; int nwiresbad = 0; for(int j=0; j<3; j++){ - int nb = std::max(0,blip->clusters[j].NWiresBad); - int nn = std::max(0,blip->clusters[j].NWiresNoisy); + int nb = std::max(0,blip.clusters[j].NWiresBad); + int nn = std::max(0,blip.clusters[j].NWiresNoisy); nwiresbad += (nb + nn); - nwirestot = blip->clusters[j].NWires; + nwirestot = blip.clusters[j].NWires; } float badwirefrac = float(nwiresbad)/nwirestot; @@ -195,52 +247,73 @@ void BlipAnalysis::addTheBlip(art::Ptr blip){ size = truncate(size,0.01); energyTrue = truncate(energyTrue,0.001); - - bool isTouchTrk = (blip->TouchTrkID >= 0 ); - _blip_id .push_back(blip->ID); - _blip_tpc .push_back(blip->TPC); - _blip_nplanes .push_back(blip->NPlanes); - _blip_proxtrkdist .push_back(blip->ProxTrkDist); - _blip_proxtrkid .push_back(blip->ProxTrkID); + bool isTouchTrk = (blip.TouchTrkID >= 0 ); + _blip_id .push_back(blip.ID); + _blip_tpc .push_back(blip.TPC); + _blip_nplanes .push_back(blip.NPlanes); + _blip_proxtrkdist .push_back(blip.ProxTrkDist); + _blip_proxtrkid .push_back(blip.ProxTrkID); _blip_touchtrk .push_back(isTouchTrk); - _blip_touchtrkid .push_back(blip->TouchTrkID); - _blip_bydeadwire .push_back(blip->clusters[2].DeadWireSep==0); + _blip_touchtrkid .push_back(blip.TouchTrkID); + _blip_bydeadwire .push_back(blip.clusters[2].DeadWireSep==0); _blip_badwirefrac .push_back(badwirefrac); _blip_x .push_back(x); _blip_y .push_back(y); _blip_z .push_back(z); - _blip_sigmayz .push_back(blip->SigmaYZ); - _blip_dx .push_back(blip->dX); - _blip_dw .push_back(blip->dYZ); + _blip_sigmayz .push_back(blip.SigmaYZ); + _blip_dx .push_back(blip.dX); + _blip_dw .push_back(blip.dYZ); _blip_size .push_back(size); _blip_charge .push_back(charge); _blip_energy .push_back(energy); - _blip_true_g4id .push_back(blip->truth.LeadG4ID); - _blip_true_pdg .push_back(blip->truth.LeadG4PDG); + _blip_true_g4id .push_back(blip.truth.LeadG4ID); + _blip_true_pdg .push_back(blip.truth.LeadG4PDG); _blip_true_energy .push_back(energyTrue); - _blip_pl0_nwires .push_back(blip->clusters[0].NWires); - _blip_pl1_nwires .push_back(blip->clusters[1].NWires); - _blip_pl2_nwires .push_back(blip->clusters[2].NWires); - _blip_pl0_nwiresbad .push_back(blip->clusters[0].NWiresBad); - _blip_pl1_nwiresbad .push_back(blip->clusters[1].NWiresBad); - _blip_pl2_nwiresbad .push_back(blip->clusters[2].NWiresBad); - _blip_pl0_bydeadwire .push_back((blip->clusters[0].DeadWireSep==0)); - _blip_pl1_bydeadwire .push_back((blip->clusters[1].DeadWireSep==0)); - _blip_pl2_bydeadwire .push_back((blip->clusters[2].DeadWireSep==0)); - _blip_pl0_centerwire .push_back(blip->clusters[0].CenterWire); - _blip_pl1_centerwire .push_back(blip->clusters[1].CenterWire); - _blip_pl2_centerwire .push_back(blip->clusters[2].CenterWire); + _blip_pl0_nwires .push_back(blip.clusters[0].NWires); + _blip_pl1_nwires .push_back(blip.clusters[1].NWires); + _blip_pl2_nwires .push_back(blip.clusters[2].NWires); + _blip_pl0_nwiresbad .push_back(blip.clusters[0].NWiresBad); + _blip_pl1_nwiresbad .push_back(blip.clusters[1].NWiresBad); + _blip_pl2_nwiresbad .push_back(blip.clusters[2].NWiresBad); + _blip_pl0_bydeadwire .push_back((blip.clusters[0].DeadWireSep==0)); + _blip_pl1_bydeadwire .push_back((blip.clusters[1].DeadWireSep==0)); + _blip_pl2_bydeadwire .push_back((blip.clusters[2].DeadWireSep==0)); + _blip_pl0_centerwire .push_back(blip.clusters[0].CenterWire); + _blip_pl1_centerwire .push_back(blip.clusters[1].CenterWire); + _blip_pl2_centerwire .push_back(blip.clusters[2].CenterWire); + if( fRerunBlipReco ) { + _blip_trkid .push_back(fBlipAlg->map_blip_trkID[blip.ID]); + _blip_trkidfrac .push_back(fBlipAlg->map_blip_trkIDfrac[blip.ID]); + _blip_rnn_dir_isValid .push_back(_blipdir_isValid); + _blip_rnn_dir_x .push_back(_blipdir.X()); + _blip_rnn_dir_y .push_back(_blipdir.Y()); + _blip_rnn_dir_z .push_back(_blipdir.Z()); + _blip_primary_ancestor_pdg .push_back(fBlipAlg->map_blip_primaryPDG[blip.ID]); + _blip_primary_ancestor_g4id.push_back(fBlipAlg->map_blip_primaryG4ID[blip.ID]); + _blip_true_category.push_back(fBlipAlg->map_blip_ncategory[blip.ID]); + TVector3 tdir(-9,-9,-9); + int g4idx = blip.truth.LeadG4Index; + if( g4idx >= 0 ) { + auto& p = fBlipAlg->pinfo[g4idx]; + tdir.SetXYZ(p.Px,p.Py,p.Pz); + tdir.SetMag(1); + } + _blip_true_dir_x.push_back(tdir.X()); + _blip_true_dir_y.push_back(tdir.Y()); + _blip_true_dir_z.push_back(tdir.Z()); + } } //---------------------------------------------------------------------------- -void BlipAnalysis::analyzeEvent(art::Event const &e, bool fData) +void BlipAnalysis::analyzeEvent(art::Event const &evt, bool fData) { - _evt = e.event(); - _sub = e.subRun(); - _run = e.run(); + _evt = evt.event(); + _sub = evt.subRun(); + _run = evt.run(); std::cout << "[BlipAnalysis::analyzeEvent] Run: " << _run << ", SubRun: " << _sub << ", Event: " << _evt << std::endl; - + + //================================================ // If a non-zero radius was configured, or we are // only saving neutrino events, then exit out of this @@ -248,24 +321,53 @@ void BlipAnalysis::analyzeEvent(art::Event const &e, bool fData) //================================================ if( fSaveOnlyNuEvts ) return; if( fNuBlipRadius > 0 ) return; - + //========================================== // ... otherwise, save all blips in the event //========================================== std::cout <<"********** BlipAnalysis::analyzeEvent **********\n"; - // reset branch vectors + resetVariables(); - // loop over blips saved to the event - art::Handle< std::vector > blipHandle; - std::vector > bliplist; - if (e.getByLabel(fBlipProducer,blipHandle)) - art::fill_ptr_vector(bliplist, blipHandle); - _nblips = bliplist.size(); - std::cout<<"Saving all "<<_nblips<<" blips in the event\n"; - for(size_t i=0; i > blipHandle; + std::vector > bliplist; + if (evt.getByLabel(fBlipProducer,blipHandle)) + art::fill_ptr_vector(bliplist, blipHandle); + _nblips = bliplist.size(); + std::cout<<"Saving all "<<_nblips<<" blips in the event\n"; + for(size_t i=0; iRunBlipReco(evt); + _elifetime = fBlipAlg->kLifetime; + _nblips = fBlipAlg->blips.size(); + for(int i=0; i<_nblips; i++){ + auto& blp = fBlipAlg->blips[i]; + _blipdir_isValid = false; + _blipdir.SetXYZ(-9,-9,-9); + if( fEnableRNN && _RNNModelLoaded ) { + std::vector dir = fBlipRNN->predict(blp,fBlipAlg->hitinfo); + if(dir.size()==3){ + _blipdir.SetXYZ(dir[0],dir[1],dir[2]); + _blipdir_isValid = ( _blipdir.Mag() < 1.01 ) ? true : false; + } + } + addTheBlip(blp); + } + } + std::cout <<"************************************************\n"; @@ -319,7 +421,7 @@ void BlipAnalysis::analyzeSlice(art::Event const &e, std::vector //TVector3 loc; TVector3 loc = (fSaveSCECorrLoc) ? blip->PositionSCE : blip->Position; if( (loc-nuvtx).Mag() > fNuBlipRadius ) continue; - addTheBlip(blip); + addTheBlip(*blip); _nblips++; } std::cout @@ -333,7 +435,8 @@ void BlipAnalysis::analyzeSlice(art::Event const &e, std::vector //---------------------------------------------------------------------------- void BlipAnalysis::setBranches(TTree *_tree) { - _tree->Branch("nblips_saved", &_nblips, "nblips_saved/I"); + _tree->Branch("elifetime", &_elifetime, "elifetime/F"); + _tree->Branch("nblips_saved", &_nblips, "nblips_saved/I"); //_tree->Branch("blip_id", "std::vector< int >", &_blip_id); _tree->Branch("blip_x", "std::vector< float >", &_blip_x); _tree->Branch("blip_y", "std::vector< float >", &_blip_y); @@ -369,6 +472,21 @@ void BlipAnalysis::setBranches(TTree *_tree) if( !fLiteMode ){ _tree->Branch("blip_true_energy", "std::vector< float >", &_blip_true_energy); } + if( fRerunBlipReco ) { + _tree->Branch("blip_trkid", "std::vector< int >", &_blip_trkid); + _tree->Branch("blip_trkidfrac", "std::vector< float >", &_blip_trkidfrac); + _tree->Branch("blip_rnn_dir_isValid","std::vector< bool >", &_blip_rnn_dir_isValid); + _tree->Branch("blip_rnn_dir_x","std::vector< float >", &_blip_rnn_dir_x); + _tree->Branch("blip_rnn_dir_y","std::vector< float >", &_blip_rnn_dir_y); + _tree->Branch("blip_rnn_dir_z","std::vector< float >", &_blip_rnn_dir_z); + _tree->Branch("blip_primary_ancestor_pdg","std::vector< int >", &_blip_primary_ancestor_pdg); + _tree->Branch("blip_primary_ancestor_g4id","std::vector< int >", &_blip_primary_ancestor_g4id); + _tree->Branch("blip_true_category","std::vector< int >", &_blip_true_category); + //_tree->Branch("blip_true_dir_x","std::vector< float >", &_blip_true_dir_x); + //_tree->Branch("blip_true_dir_y","std::vector< float >", &_blip_true_dir_y); + //_tree->Branch("blip_true_dir_z","std::vector< float >", &_blip_true_dir_z); + + } } @@ -382,7 +500,8 @@ void BlipAnalysis::resetTTree(TTree *_tree) //---------------------------------------------------------------------------- void BlipAnalysis::resetVariables() { - _nblips = 0; + _nblips = 0; + _elifetime = 0; _blip_id .clear(); _blip_tpc .clear(); _blip_nplanes .clear(); @@ -416,6 +535,18 @@ void BlipAnalysis::resetVariables() _blip_pl0_centerwire.clear(); _blip_pl1_centerwire.clear(); _blip_pl2_centerwire.clear(); + _blip_trkid .clear(); + _blip_trkidfrac .clear(); + _blip_rnn_dir_isValid.clear(); + _blip_rnn_dir_x .clear(); + _blip_rnn_dir_y .clear(); + _blip_rnn_dir_z .clear(); + _blip_primary_ancestor_pdg.clear(); + _blip_primary_ancestor_g4id.clear(); + _blip_true_category.clear(); + _blip_true_dir_x .clear(); + _blip_true_dir_y .clear(); + _blip_true_dir_z .clear(); } diff --git a/ubana/searchingfornues/Selection/AnalysisTools/CMakeLists.txt b/ubana/searchingfornues/Selection/AnalysisTools/CMakeLists.txt index 63e23b85e..26f866d28 100644 --- a/ubana/searchingfornues/Selection/AnalysisTools/CMakeLists.txt +++ b/ubana/searchingfornues/Selection/AnalysisTools/CMakeLists.txt @@ -352,6 +352,7 @@ cet_build_plugin( nusimdata::SimulationBase ROOT::EG ubreco::BlipReco_Alg + ubreco::BlipReco_RNN ) cet_build_plugin( diff --git a/ubana/searchingfornues/Selection/job/neutrinoselection.fcl b/ubana/searchingfornues/Selection/job/neutrinoselection.fcl index 48575f8a8..5575000e3 100644 --- a/ubana/searchingfornues/Selection/job/neutrinoselection.fcl +++ b/ubana/searchingfornues/Selection/job/neutrinoselection.fcl @@ -3,6 +3,7 @@ #include "proximityclustering.fcl" #include "eventweight_microboone_genie_knobs.fcl" + #include "microboone_blipreco.fcl" FidVol: { @@ -90,7 +91,15 @@ BlipAnalysisTool: { NuBlipRadius: -1 # Only save blips within this dist of nu vtx [cm] # ^^ set <= 0 to save ALL blips in neutrino events LiteMode: false # save bare minimum (x,y,z,size,energy,trueG4ID) + RerunBlipReco: true # re-run blip reconstruction instead of reading saved reco2-blips + BlipAlg: @local::microboone_blipalg + EnableRNN: true # enable RNN direction predictor for proton blips + # (^^ REQUIRES re-run of blipreco) } +BlipAnalysisTool.BlipAlg.CalAreaConstants: [4.31e-3, 4.02e-3, 4.10e-3] +BlipAnalysisTool.BlipAlg.DebugMode: false +BlipAnalysisTool.BlipRNN.ModelFilename: "model_blipRNN.pt" + BlipAnalysisToolMCC9: { tool_type: "BlipAnalysisMCC9"