diff --git a/sbncode/CAFMaker/CAFMakerParams.h b/sbncode/CAFMaker/CAFMakerParams.h index 271e54353..92188f8af 100644 --- a/sbncode/CAFMaker/CAFMakerParams.h +++ b/sbncode/CAFMaker/CAFMakerParams.h @@ -204,6 +204,12 @@ namespace caf "pandoraShowerRazzle" }; + Atom PFPRazzledLabel { + Name("PFPRazzledLabel"), + Comment("Base label of pfp mva particle-id producer."), + "pandoraRazzled" + }; + Atom RecoShowerSelectionLabel { Name("RecoShowerSelectionLabel"), Comment("Base label of shower selection vars producer."), diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 7f58a46dc..439eb9857 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -1686,6 +1686,10 @@ void CAFMaker::produce(art::Event& evt) noexcept { FindManyPStrict(slcShowers, evt, fParams.ShowerRazzleLabel() + slice_tag_suff); + art::FindManyP fmPFPRazzled = + FindManyPStrict(fmPFPart, evt, + fParams.PFPRazzledLabel() + slice_tag_suff); + art::FindOneP foCNNScores = FindOnePStrict(fmPFPart, evt, fParams.CNNScoreLabel() + slice_tag_suff); @@ -1916,6 +1920,10 @@ void CAFMaker::produce(art::Event& evt) noexcept { } } + if (fmPFPRazzled.isValid() && fmPFPRazzled.at(iPart).size()==1) { + FillPFPRazzled(fmPFPRazzled.at(iPart).front(), pfp); + } + // fill all the stuff SRTrack& trk = pfp.trk; FillTrackVars(*thisTrack[0], producer, trk); diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index 548dc4c49..095cf8028 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/sbncode/CAFMaker/FillReco.cxx @@ -996,6 +996,22 @@ namespace caf //...................................................................... + void FillPFPRazzled(const art::Ptr razzled, + caf::SRPFP& srpfp, + bool allowEmpty) + { + srpfp.razzled.electronScore = razzled->mvaScoreMap.at(11); + srpfp.razzled.muonScore = razzled->mvaScoreMap.at(13); + srpfp.razzled.photonScore = razzled->mvaScoreMap.at(22); + srpfp.razzled.pionScore = razzled->mvaScoreMap.at(211); + srpfp.razzled.protonScore = razzled->mvaScoreMap.at(2212); + + srpfp.razzled.pdg = razzled->BestPDG(); + srpfp.razzled.bestScore = razzled->BestScore(); + } + + //...................................................................... + void SetNuMuCCPrimary(std::vector &recs, std::vector &srneutrinos) { // // set is_primary to true by default diff --git a/sbncode/CAFMaker/FillReco.h b/sbncode/CAFMaker/FillReco.h index a90de4ada..fb40df775 100644 --- a/sbncode/CAFMaker/FillReco.h +++ b/sbncode/CAFMaker/FillReco.h @@ -231,6 +231,10 @@ namespace caf void FillTPCPMTBarycenterMatch(const sbn::TPCPMTBarycenterMatch *matchInfo, caf::SRSlice& slice); + void FillPFPRazzled(const art::Ptr razzled, + caf::SRPFP& srpfp, + bool allowEmpty = false); + template void CopyPropertyIfSet( const std::map& props, const std::string& search, U& value ); } diff --git a/sbncode/PID/CMakeLists.txt b/sbncode/PID/CMakeLists.txt index 1ad88d077..1d5da7ad7 100644 --- a/sbncode/PID/CMakeLists.txt +++ b/sbncode/PID/CMakeLists.txt @@ -53,6 +53,19 @@ cet_build_plugin( Razzle art::module caf_RecoUtils ) +cet_build_plugin( Razzled art::module + LIBRARIES + ROOT::TMVA + canvas::canvas + art_root_io::TFileService_service + art::Persistency_Provenance + larsim::Utils + larsim::MCCheater_BackTrackerService_service + larsim::MCCheater_ParticleInventoryService_service + lardata::Utilities + sbnobj::Common_Reco +) + install_headers() install_fhicl() install_source() diff --git a/sbncode/PID/Razzle_module.cc b/sbncode/PID/Razzle_module.cc index 70c469ecd..498e57eff 100644 --- a/sbncode/PID/Razzle_module.cc +++ b/sbncode/PID/Razzle_module.cc @@ -479,12 +479,12 @@ void Razzle::FillShowerMetrics(const recob::Shower& shower, const std::vector 1st induction > 2nd induction - bestPlane = shower.best_plane(); - - if (bestPlane < 0 || bestPlane > 3) - throw cet::exception("Razzle") << "Best plane: " << bestPlane; + if(showerPlaneHits[2] >= showerPlaneHits[1] && showerPlaneHits[2] >= showerPlaneHits[0]) + bestPlane = 2; + else if(showerPlaneHits[0] >= showerPlaneHits[1]) + bestPlane = 0; + else + bestPlane = 1; bestdEdx = shower.dEdx()[bestPlane]; bestdEdx = std::min(bestdEdx, 20.f); diff --git a/sbncode/PID/Razzled_module.cc b/sbncode/PID/Razzled_module.cc new file mode 100644 index 000000000..e7bafd699 --- /dev/null +++ b/sbncode/PID/Razzled_module.cc @@ -0,0 +1,958 @@ +//////////////////////////////////////////////////////////////////////// +// Class: Razzled +// Plugin Type: producer (art v3_05_01) +// File: Razzled_module.cc +// +// Created by Henry Lay, May 2023 to combine the work done by Edward +// Tyley in the Razzle and Dazzle modules into a single tool for +// PFP PID regardless of track or shower decision. +// +// Module that attempts to classify each shower as either: +// - Electron (11) +// - Muon (13) +// - Photon (22) +// - Pion (211) +// - Proton (2212) +// +// If a PFP is not classified it will be assigned a pdg of -5 +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/FindOneP.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "art_root_io/TFileService.h" + +#include "lardata/Utilities/AssociationUtil.h" +#include "lardata/DetectorInfoServices/DetectorClocksService.h" + +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/PFParticle.h" +#include "lardataobj/RecoBase/PFParticleMetadata.h" +#include "lardataobj/RecoBase/Cluster.h" +#include "lardataobj/RecoBase/Track.h" +#include "lardataobj/RecoBase/Shower.h" +#include "lardataobj/RecoBase/Vertex.h" +#include "lardataobj/RecoBase/MCSFitResult.h" +#include "lardataobj/AnalysisBase/Calorimetry.h" +#include "lardataobj/AnalysisBase/ParticleID.h" + +#include "larcore/CoreUtils/ServiceUtil.h" +#include "larcore/Geometry/Geometry.h" + +#include "larsim/MCCheater/BackTrackerService.h" +#include "larsim/MCCheater/ParticleInventoryService.h" +#include "larsim/Utils/TruthMatchUtils.h" + +#include "sbnobj/Common/Reco/MVAPID.h" +#include "sbnobj/Common/Reco/RangeP.h" +#include "sbnobj/Common/Reco/ScatterClosestApproach.h" +#include "sbnobj/Common/Reco/StoppingChi2Fit.h" + +#include "TTree.h" +#include "TMVA/Reader.h" + +#include + +namespace sbn { + class Razzled : public art::EDProducer { + public: + explicit Razzled(fhicl::ParameterSet const& p); + + Razzled(Razzled const&) = delete; + Razzled(Razzled&&) = delete; + Razzled& operator=(Razzled const&) = delete; + Razzled& operator=(Razzled&&) = delete; + + void produce(art::Event& e) override; + void beginJob() override; + + private: + art::ServiceHandle tfs; + + art::InputTag fPFPLabel, fClusterLabel, fTrackLabel, fShowerLabel, fCaloLabel, + fMCSLabel, fChi2Label, fRangeLabel, fClosestApproachLabel, fStoppingChi2Label, + fDazzleLabel, fRazzleLabel; + + const float fMinTrackLength, fMinShowerEnergy; + const bool fMakeTree, fRunMVA, fSaveFullCalo; + const std::string fMethodName, fWeightFile; + const float fXMin, fXMax, fYMin, fYMax, fZMin, fZMax; + + float pfp_numDaughters; // PFP's number of daughters in particle hierarchy + float pfp_maxDaughterHits; // Max number of hits from any PFP daughters + float pfp_trackScore; + float pfp_chargeEndFrac; + float pfp_chargeFracSpread; + float pfp_linearFitDiff; + float pfp_linearFitLength; + float pfp_linearFitGapLength; + float pfp_linearFitRMS; + float pfp_openAngleDiff; + float pfp_secondaryPCARatio; + float pfp_tertiaryPCARatio; + float pfp_vertexDist; + + float trk_length; // The length of the track [cm] + float trk_chi2PIDMuon; // Chi2 PID score for muon hypothesis + float trk_chi2PIDProton; // Chi2 PID score for proton hypothesis + float trk_chi2PIDMuonPionDiff; // The muon and pion scores are very correlated, instead take the relative difference + float trk_mcsScatterMean; // Mean scattering angle from MCS + float trk_mcsScatterMaxRatio; // Scattering angles used in MCS calculation [mrad] + float trk_momDiff; // Relative momentum agreement between range and MCS + float trk_meanDCA; // Distance of closest approach to interpolated track [cm] + float trk_stoppingdEdxChi2Ratio; // Ratio of exp/pol0 chi2 fits to end of track + float trk_chi2Pol0dEdxFit; // Fitted pol0 to find the dE/dx of the track [MeV/cm] + + float shw_bestdEdx; // The dE/dx at the start of the shower (in the best plane) [MeV/cm] + float shw_convGap; // The gap between the shower start and parent vertex [cm] + float shw_openAngle; // Opening Angle of the shower, defined as atan(width/length) [deg] + float shw_modHitDensity; // The hit density corrected for the pitch (hits/wire) + float shw_sqrtEnergyDensity; // Sqrt of the energy divided by the length [cm^-1] + + + float electronScore, muonScore, photonScore, pionScore, protonScore, bestScore; + int bestPDG; + + TMVA::Reader* reader; + TTree* pfpTree; + + int truePDG, trueMotherPDG; + std::string trueType, trueEndProcess; + float energyComp, energyPurity, trueEndMomentum; + + bool recoPrimary, trackContained, showerContained, unambiguousSlice, goodTrack, goodShower; + int recoPDG, chi2PDG; + float trackStartX, trackStartY, trackStartZ, trackEndX, trackEndY, trackEndZ, trackChi2PIDPion, trackChi2PIDKaon, + showerStartX, showerStartY, showerStartZ, showerEndX, showerEndY, showerEndZ, + showerEnergy; + + float dazzleMuonScore, dazzlePionScore, dazzleProtonScore, dazzleOtherScore, + razzleElectronScore, razzlePhotonScore, razzleOtherScore; + int dazzlePDG, razzlePDG; + + std::vector trackdEdx, trackResRange; + + void ClearTreeValues(); + + void FillTrueParticleMetrics(const detinfo::DetectorClocksData &clockData, const std::vector> &hits); + + void FillPFPMetrics(const art::Ptr &pfp, const std::map> &pfpMap, + const art::FindManyP &pfpsToClusters, const art::FindManyP &clustersToHits, + const art::FindOneP &pfpsToMetadata); + + void FillPandoraTrackIDScoreVars(const std::map &propertiesMap); + + void FillPandoraTrackIDScoreVar(const std::map &propertiesMap, float &var, std::string name); + + void FillTrackMetrics(const art::Ptr &track); + + void FillMCSMetrics(const art::Ptr &mcs); + + void FillMomDiff(const art::Ptr &rangeP, const art::Ptr &mcs); + + void FillChi2PIDMetrics(const art::Ptr &pid); + + void FillStoppingChi2Metrics(const art::Ptr &stoppingChi2); + + void FillDazzleMetrics(const art::Ptr &dazzle); + + void FillShowerMetrics(const art::Ptr &shower, const std::vector> &hits); + + void FillShowerConversionGap(const art::Ptr &pfp, const std::map> &pfpMap, + const art::Ptr &shower, const art::FindOneP &pfpsToVertices); + + void FillRazzleMetrics(const art::Ptr &razzle); + + bool InFV(const TVector3 &pos); + + std::map> GetPFPMap(const std::vector> &pfps); + + std::string PDGString(const int pdg); + + MVAPID RunMVA(); + }; + + Razzled::Razzled(fhicl::ParameterSet const& p) + : EDProducer { p } + , fPFPLabel(p.get("PFPLabel")) + , fClusterLabel(p.get("ClusterLabel")) + , fTrackLabel(p.get("TrackLabel")) + , fShowerLabel(p.get("ShowerLabel")) + , fCaloLabel(p.get("CaloLabel")) + , fMCSLabel(p.get("MCSLabel"), std::string("muon")) + , fChi2Label(p.get("Chi2Label")) + , fRangeLabel(p.get("RangeLabel"), std::string("muon")) + , fClosestApproachLabel(p.get("ClosestApproachLabel")) + , fStoppingChi2Label(p.get("StoppingChi2Label")) + , fDazzleLabel(p.get("DazzleLabel")) + , fRazzleLabel(p.get("RazzleLabel")) + , fMinTrackLength(p.get("MinTrackLength")) + , fMinShowerEnergy(p.get("MinShowerEnergy")) + , fMakeTree(p.get("MakeTree")) + , fRunMVA(p.get("RunMVA")) + , fSaveFullCalo(p.get("SaveFullCalo")) + , fMethodName(p.get("MethodName", "")) + , fWeightFile(p.get("WeightFile", "")) + , fXMin(p.get("XMin")) + , fXMax(p.get("XMax")) + , fYMin(p.get("YMin")) + , fYMax(p.get("YMax")) + , fZMin(p.get("ZMin")) + , fZMax(p.get("ZMax")) + { + if(!fMakeTree && !fRunMVA) + throw cet::exception("Razzled") << "Configured to do nothing"; + + if(fRunMVA) + { + if(fMethodName == "" || fWeightFile == "") + throw cet::exception("Razzled") << "Trying to run MVA with inputs not set: MethodName: " << fMethodName << " and WeightFile: " << fWeightFile; + + cet::search_path searchPath("FW_SEARCH_PATH"); + std::string fWeightFileFullPath; + if(!searchPath.find_file(fWeightFile, fWeightFileFullPath)) + throw cet::exception("Razzled") << "Unable to find weight file: " << fWeightFile << " in FW_SEARCH_PATH: " << searchPath.to_string(); + + reader = new TMVA::Reader("V"); + + reader->AddVariable("pfp_numDaughters", &pfp_numDaughters); + reader->AddVariable("pfp_maxDaughterHits", &pfp_maxDaughterHits); + reader->AddVariable("pfp_trackScore", &pfp_trackScore); + + reader->AddVariable("trk_length", &trk_length); + reader->AddVariable("trk_chi2PIDMuon", &trk_chi2PIDMuon); + reader->AddVariable("trk_chi2PIDProton", &trk_chi2PIDProton); + reader->AddVariable("trk_chi2PIDMuonPionDiff", &trk_chi2PIDMuonPionDiff); + reader->AddVariable("trk_mcsScatterMean", &trk_mcsScatterMean); + reader->AddVariable("trk_mcsScatterMaxRatio", &trk_mcsScatterMaxRatio); + reader->AddVariable("trk_meanDCA", &trk_meanDCA); + reader->AddVariable("trk_stoppingdEdxChi2Ratio", &trk_stoppingdEdxChi2Ratio); + reader->AddVariable("trk_chi2Pol0dEdxFit", &trk_chi2Pol0dEdxFit); + reader->AddVariable("trk_momDiff", &trk_momDiff); + + reader->AddVariable("shw_bestdEdx", &shw_bestdEdx); + reader->AddVariable("shw_convGap", &shw_convGap); + reader->AddVariable("shw_openAngle", &shw_openAngle); + reader->AddVariable("shw_modHitDensity", &shw_modHitDensity); + reader->AddVariable("shw_sqrtEnergyDensity>2.5?2.5:shw_sqrtEnergyDensity", &shw_sqrtEnergyDensity); + + reader->BookMVA(fMethodName, fWeightFileFullPath); + } + + produces>(); + produces>(); + } + + void Razzled::beginJob() + { + if(fMakeTree) + { + pfpTree = tfs->make("pfpTree", "Tree filled per PFP with PID variables"); + + if(fRunMVA) + { + pfpTree->Branch("electronScore", &electronScore); + pfpTree->Branch("muonScore", &muonScore); + pfpTree->Branch("photonScore", &photonScore); + pfpTree->Branch("pionScore", &pionScore); + pfpTree->Branch("protonScore", &protonScore); + pfpTree->Branch("bestScore", &bestScore); + pfpTree->Branch("bestPDG", &bestPDG); + } + + pfpTree->Branch("truePDG", &truePDG); + pfpTree->Branch("trueMotherPDG", &trueMotherPDG); + pfpTree->Branch("trueType", &trueType); + pfpTree->Branch("trueEndProcess", &trueEndProcess); + pfpTree->Branch("trueEndMomentum", &trueEndMomentum); + pfpTree->Branch("energyComp", &energyComp); + pfpTree->Branch("energyPurity", &energyPurity); + + pfpTree->Branch("recoPrimary", &recoPrimary); + pfpTree->Branch("unambiguousSlice", &unambiguousSlice); + pfpTree->Branch("recoPDG", &recoPDG); + + pfpTree->Branch("trackStartX", &trackStartX); + pfpTree->Branch("trackStartY", &trackStartY); + pfpTree->Branch("trackStartZ", &trackStartZ); + pfpTree->Branch("trackEndX", &trackEndX); + pfpTree->Branch("trackEndY", &trackEndY); + pfpTree->Branch("trackEndZ", &trackEndZ); + pfpTree->Branch("trackChi2PIDPion", &trackChi2PIDPion); + pfpTree->Branch("trackChi2PIDKaon", &trackChi2PIDKaon); + pfpTree->Branch("chi2PDG", &chi2PDG); + pfpTree->Branch("dazzleMuonScore", &dazzleMuonScore); + pfpTree->Branch("dazzlePionScore", &dazzlePionScore); + pfpTree->Branch("dazzleProtonScore", &dazzleProtonScore); + pfpTree->Branch("dazzleOtherScore", &dazzleOtherScore); + pfpTree->Branch("dazzlePDG", &dazzlePDG); + pfpTree->Branch("trackContained", &trackContained); + pfpTree->Branch("goodTrack", &goodTrack); + + pfpTree->Branch("showerStartX", &showerStartX); + pfpTree->Branch("showerStartY", &showerStartY); + pfpTree->Branch("showerStartZ", &showerStartZ); + pfpTree->Branch("showerEndX", &showerEndX); + pfpTree->Branch("showerEndY", &showerEndY); + pfpTree->Branch("showerEndZ", &showerEndZ); + pfpTree->Branch("showerEnergy", &showerEnergy); + pfpTree->Branch("razzleElectronScore", &razzleElectronScore); + pfpTree->Branch("razzlePhotonScore", &razzlePhotonScore); + pfpTree->Branch("razzleOtherScore", &razzleOtherScore); + pfpTree->Branch("razzlePDG", &razzlePDG); + pfpTree->Branch("showerContained", &showerContained); + pfpTree->Branch("goodShower", &goodShower); + + pfpTree->Branch("pfp_numDaughters", &pfp_numDaughters); + pfpTree->Branch("pfp_maxDaughterHits", &pfp_maxDaughterHits); + pfpTree->Branch("pfp_trackScore", &pfp_trackScore); + pfpTree->Branch("pfp_chargeEndFrac", &pfp_chargeEndFrac); + pfpTree->Branch("pfp_chargeFracSpread", &pfp_chargeFracSpread); + pfpTree->Branch("pfp_linearFitDiff", &pfp_linearFitDiff); + pfpTree->Branch("pfp_linearFitLength", &pfp_linearFitLength); + pfpTree->Branch("pfp_linearFitGapLength", &pfp_linearFitGapLength); + pfpTree->Branch("pfp_linearFitRMS", &pfp_linearFitRMS); + pfpTree->Branch("pfp_openAngleDiff", &pfp_openAngleDiff); + pfpTree->Branch("pfp_secondaryPCARatio", &pfp_secondaryPCARatio); + pfpTree->Branch("pfp_tertiaryPCARatio", &pfp_tertiaryPCARatio); + pfpTree->Branch("pfp_vertexDist", &pfp_vertexDist); + + pfpTree->Branch("trk_length", &trk_length); + pfpTree->Branch("trk_chi2PIDMuon", &trk_chi2PIDMuon); + pfpTree->Branch("trk_chi2PIDProton", &trk_chi2PIDProton); + pfpTree->Branch("trk_chi2PIDMuonPionDiff", &trk_chi2PIDMuonPionDiff); + pfpTree->Branch("trk_mcsScatterMean", &trk_mcsScatterMean); + pfpTree->Branch("trk_mcsScatterMaxRatio", &trk_mcsScatterMaxRatio); + pfpTree->Branch("trk_meanDCA", &trk_meanDCA); + pfpTree->Branch("trk_stoppingdEdxChi2Ratio", &trk_stoppingdEdxChi2Ratio); + pfpTree->Branch("trk_chi2Pol0dEdxFit", &trk_chi2Pol0dEdxFit); + pfpTree->Branch("trk_momDiff", &trk_momDiff); + + pfpTree->Branch("shw_bestdEdx", &shw_bestdEdx); + pfpTree->Branch("shw_convGap", &shw_convGap); + pfpTree->Branch("shw_openAngle", &shw_openAngle); + pfpTree->Branch("shw_modHitDensity", &shw_modHitDensity); + pfpTree->Branch("shw_sqrtEnergyDensity", &shw_sqrtEnergyDensity); + + if(fSaveFullCalo) + { + pfpTree->Branch("trackdEdx", &trackdEdx); + pfpTree->Branch("trackResRange", &trackResRange); + } + } + } + + void Razzled::produce(art::Event& e) + { + auto const clockData(art::ServiceHandle()->DataFor(e)); + + auto const pfpHandle(e.getValidHandle>(fPFPLabel)); + auto const clusterHandle(e.getValidHandle>(fClusterLabel)); + auto const trackHandle(e.getValidHandle>(fTrackLabel)); + auto const showerHandle(e.getValidHandle>(fShowerLabel)); + + std::vector> pfps; + art::fill_ptr_vector(pfps, pfpHandle); + + art::FindOneP pfpsToVertices(pfpHandle, e, fPFPLabel); + art::FindOneP pfpsToMetadata(pfpHandle, e, fPFPLabel); + art::FindManyP pfpsToClusters(pfpHandle, e, fPFPLabel); + art::FindManyP clustersToHits(clusterHandle, e, fClusterLabel); + + art::FindOneP pfpsToTracks(pfpHandle, e, fTrackLabel); + art::FindManyP tracksToCalos(trackHandle, e, fCaloLabel); + art::FindOneP tracksToMCSs(trackHandle, e, fMCSLabel); + art::FindManyP tracksToChi2s(trackHandle, e, fChi2Label); + art::FindOneP tracksToRangePs(trackHandle, e, fRangeLabel); + art::FindOneP tracksToClosestApproaches(trackHandle, e, fClosestApproachLabel); + art::FindOneP tracksToStoppingChi2s(trackHandle, e, fStoppingChi2Label); + art::FindOneP tracksToDazzles(trackHandle, e, fDazzleLabel); + + art::FindOneP pfpsToShowers(pfpHandle, e, fShowerLabel); + art::FindManyP showersToHits(showerHandle, e, fShowerLabel); + art::FindOneP showersToRazzles(showerHandle, e, fRazzleLabel); + + auto mvaPIDVec = std::make_unique>(); + auto pfpAssns = std::make_unique>(); + + const std::map> pfpMap = this->GetPFPMap(pfps); + + for(auto const& pfp : pfps) + { + this->ClearTreeValues(); + + const art::Ptr track = pfpsToTracks.at(pfp.key()); + const art::Ptr shower = pfpsToShowers.at(pfp.key()); + + if(track.isNull() && shower.isNull()) + continue; + + this->FillPFPMetrics(pfp, pfpMap, pfpsToClusters, clustersToHits, pfpsToMetadata); + + if(track.isNonnull()) + { + goodTrack = true; + + if(track->Length() < fMinTrackLength) + continue; + + this->FillTrackMetrics(track); + + const std::vector> calos = tracksToCalos.at(track.key()); + if(calos.size() != 3) + continue; + + const unsigned int maxHits = std::max({ calos[0]->dEdx().size(), calos[1]->dEdx().size(), calos[2]->dEdx().size()}); + const int bestPlane = (calos[2]->dEdx().size() == maxHits) ? 2 : (calos[0]->dEdx().size() == maxHits) ? 0 : (calos[1]->dEdx().size() == maxHits) ? 1 : -1; + + if(bestPlane < 0 || bestPlane > 3) + throw cet::exception("Razzled") << "Best plane: " << bestPlane; + + if(fSaveFullCalo) + { + trackdEdx = calos[bestPlane]->dEdx(); + trackResRange = calos[bestPlane]->ResidualRange(); + } + + const std::vector> chi2s = tracksToChi2s.at(track.key()); + if(chi2s.size() == 3) + this->FillChi2PIDMetrics(chi2s[bestPlane]); + + const art::Ptr mcs = tracksToMCSs.at(track.key()); + if(mcs.isNonnull()) + this->FillMCSMetrics(mcs); + + const art::Ptr rangeP = tracksToRangePs.at(track.key()); + if(rangeP.isNonnull() && mcs.isNonnull()) + this->FillMomDiff(rangeP, mcs); + + const art::Ptr closestApproach = tracksToClosestApproaches.at(track.key()); + if(closestApproach.isNonnull()) + trk_meanDCA = closestApproach->mean; + + const art::Ptr stoppingChi2 = tracksToStoppingChi2s.at(track.key()); + if(stoppingChi2.isNonnull()) + this->FillStoppingChi2Metrics(stoppingChi2); + + const art::Ptr dazzle = tracksToDazzles.at(track.key()); + if(dazzle.isNonnull()) + this->FillDazzleMetrics(dazzle); + } + + if(shower.isNonnull()) + { + goodShower = true; + + if(shower->best_plane() < 0 || shower->Energy().at(shower->best_plane()) < fMinShowerEnergy) + continue; + + showerEnergy = shower->Energy().at(shower->best_plane()); + + const std::vector> showerHits = showersToHits.at(shower.key()); + + this->FillShowerMetrics(shower, showerHits); + this->FillShowerConversionGap(pfp, pfpMap, shower, pfpsToVertices); + + const art::Ptr razzle = showersToRazzles.at(shower.key()); + if(razzle.isNonnull()) + this->FillRazzleMetrics(razzle); + } + + if(fRunMVA) + { + MVAPID mvaPID = this->RunMVA(); + mvaPIDVec->push_back(mvaPID); + util::CreateAssn(*this, e, *mvaPIDVec, pfp, *pfpAssns); + } + + if(fMakeTree) + { + std::vector> hits; + + const std::vector> clusters = pfpsToClusters.at(pfp.key()); + for(auto const& cluster : clusters) + { + const std::vector> clusterHits = clustersToHits.at(cluster.key()); + hits.insert(hits.end(), clusterHits.begin(), clusterHits.end()); + } + + this->FillTrueParticleMetrics(clockData, hits); + pfpTree->Fill(); + } + } + + e.put(std::move(mvaPIDVec)); + e.put(std::move(pfpAssns)); + } + + void Razzled::ClearTreeValues() + { + muonScore = -5.f; electronScore = -5.f; photonScore = -5.f; pionScore = -5.f; + protonScore = -5.f; bestScore = -5.f; + bestPDG = -5; + + truePDG = -5; trueMotherPDG = -5; trueType = ""; trueEndProcess = ""; trueEndMomentum = -5.f; + energyComp = -5.f; energyPurity = -5.f; + + recoPrimary = false; unambiguousSlice = false; recoPDG = -5; + + trackStartX = -999.f; trackStartY = -999.f; trackStartZ = -999.f; + trackEndX = -999.f; trackEndY = -999.f; trackEndZ = -999.f; + trackChi2PIDPion = -999.f; trackChi2PIDKaon = -999.f; + chi2PDG = -1; + trackContained = false; goodTrack = false; + + dazzleMuonScore = -2.f; dazzlePionScore = -2.f; dazzleProtonScore = -2.f; + dazzleOtherScore = -2.f; + dazzlePDG = -1; + + showerStartX = -999.f; showerStartY = -999.f; showerStartZ = -999.f; + showerEndX = -999.f; showerEndY = -999.f; showerEndZ = -999.f; + showerEnergy = -999.f; + showerContained = false; goodShower = false; + + razzleElectronScore = -2.f; razzlePhotonScore = -2.f; razzleOtherScore = -2.f; + razzlePDG = -1; + + pfp_numDaughters = -5.f; pfp_maxDaughterHits = -50.f; pfp_trackScore = -1.f; + pfp_chargeEndFrac = -1.f; pfp_chargeFracSpread = -1.f; pfp_linearFitDiff = -1.f; + pfp_linearFitLength = -50.f; pfp_linearFitGapLength = -1.f; pfp_linearFitRMS = -1.f; + pfp_openAngleDiff = -10.f; pfp_secondaryPCARatio = -1.f; pfp_tertiaryPCARatio = -1.f; + pfp_vertexDist = -50.f; + + trk_length = -100.f; trk_chi2PIDMuon = -10.f; trk_chi2PIDProton = -10.f; + trk_chi2PIDMuonPionDiff = -100.f; trk_mcsScatterMean = -100.f; trk_mcsScatterMaxRatio = -1.f; + trk_meanDCA = -5.f; trk_stoppingdEdxChi2Ratio = -5.f; trk_chi2Pol0dEdxFit = -5.f; + trk_momDiff = -10.f; + + shw_bestdEdx = -5.f; shw_convGap = -5.f; shw_openAngle = -10.f; + shw_modHitDensity = -5.f; shw_sqrtEnergyDensity = -5.f; + + trackdEdx.clear(); trackResRange.clear(); + } + + void Razzled::FillTrueParticleMetrics(const detinfo::DetectorClocksData &clockData, const std::vector> &hits) + { + art::ServiceHandle bt_serv; + art::ServiceHandle particleInv; + + const int bestMatch = TruthMatchUtils::TrueParticleIDFromTotalTrueEnergy(clockData, hits, true) ; + + if(!TruthMatchUtils::Valid(bestMatch)) + return; + + float totalHitEnergy = 0.f, totalTrueHitEnergy = 0.f; + + for(auto const& hit : hits) + { + const std::vector ides = bt_serv->HitToTrackIDEs(clockData, hit); + + totalHitEnergy = std::accumulate(ides.cbegin(), ides.cend(), totalHitEnergy, + [](float sum, auto const& ide) { return sum + ide.energy; }); + totalTrueHitEnergy = std::accumulate(ides.cbegin(), ides.cend(), totalTrueHitEnergy, + [bestMatch](float sum, auto const& ide) { return (std::abs(ide.trackID) == bestMatch) ? sum + ide.energy : sum; }); + } + + const std::vector allIDEs = bt_serv->TrackIdToSimIDEs_Ps(bestMatch); + + float totalTrueEnergy = std::accumulate(allIDEs.cbegin(), allIDEs.cend(), 0.f, + [](float sum, auto const& ide) { return sum + ide->energy; }); + + energyComp = totalTrueHitEnergy / totalTrueEnergy; + energyPurity = totalTrueHitEnergy / totalHitEnergy; + + const simb::MCParticle* trueParticle = particleInv->TrackIdToParticle_P(bestMatch); + + if(trueParticle == NULL) + return; + + truePDG = trueParticle->PdgCode(); + trueType = this->PDGString(truePDG); + trueEndProcess = trueParticle->EndProcess(); + trueEndMomentum = trueParticle->EndMomentum().Vect().Mag(); + + if(trueParticle->Process() == "primary") + return; + + const simb::MCParticle* trueMother = particleInv->TrackIdToParticle_P(trueParticle->Mother()); + + if(trueMother == NULL) + trueMotherPDG = 0; + else + trueMotherPDG = trueMother->PdgCode(); + } + + void Razzled::FillPFPMetrics(const art::Ptr &pfp, const std::map> &pfpMap, + const art::FindManyP &pfpsToClusters, const art::FindManyP &clustersToHits, + const art::FindOneP &pfpsToMetadata) + { + recoPDG = pfp->PdgCode(); + pfp_numDaughters = pfp->Daughters().size(); + + for(auto const& daughterID : pfp->Daughters()) + { + auto const &daughterIter = pfpMap.find(daughterID); + if(daughterIter == pfpMap.end()) + continue; + + int hits = 0; + + const std::vector> clusters = pfpsToClusters.at(daughterIter->second.key()); + for(auto const& cluster : clusters) + { + const std::vector> clusterHits = clustersToHits.at(cluster.key()); + hits += clusterHits.size(); + } + + if(hits > pfp_maxDaughterHits) + pfp_maxDaughterHits = hits; + } + + const art::Ptr metadata(pfpsToMetadata.at(pfp.key())); + + if(metadata.isNull()) + return; + + const std::map propertiesMap = metadata->GetPropertiesMap(); + this->FillPandoraTrackIDScoreVars(propertiesMap); + + if(pfp->IsPrimary()) + { + recoPrimary = true; + unambiguousSlice = pfp->PdgCode() == 13 || pfp->PdgCode() == 11; + return; + } + + auto const& parentIter = pfpMap.find(pfp->Parent()); + + if(parentIter == pfpMap.end()) + return; + + art::Ptr parent = parentIter->second; + + recoPrimary = parent->IsPrimary(); + + while(!parent->IsPrimary()) + { + auto const& parentIter2 = pfpMap.find(parent->Parent()); + + if(parentIter2 == pfpMap.end()) + return; + + parent = parentIter2->second; + } + + unambiguousSlice = parent->PdgCode() == 13 || parent->PdgCode() == 11; + } + + void Razzled::FillPandoraTrackIDScoreVars(const std::map &propertiesMap) + { + this->FillPandoraTrackIDScoreVar(propertiesMap, pfp_trackScore, "TrackScore"); + this->FillPandoraTrackIDScoreVar(propertiesMap, pfp_chargeEndFrac, "LArThreeDChargeFeatureTool_EndFraction"); + this->FillPandoraTrackIDScoreVar(propertiesMap, pfp_chargeFracSpread, "LArThreeDChargeFeatureTool_FractionalSpread"); + this->FillPandoraTrackIDScoreVar(propertiesMap, pfp_linearFitDiff, "LArThreeDLinearFitFeatureTool_DiffStraightLineMean"); + this->FillPandoraTrackIDScoreVar(propertiesMap, pfp_linearFitLength, "LArThreeDLinearFitFeatureTool_Length"); + this->FillPandoraTrackIDScoreVar(propertiesMap, pfp_linearFitGapLength, "LArThreeDLinearFitFeatureTool_MaxFitGapLength"); + this->FillPandoraTrackIDScoreVar(propertiesMap, pfp_linearFitRMS, "LArThreeDLinearFitFeatureTool_SlidingLinearFitRMS"); + this->FillPandoraTrackIDScoreVar(propertiesMap, pfp_openAngleDiff, "LArThreeDOpeningAngleFeatureTool_AngleDiff"); + pfp_openAngleDiff *= TMath::RadToDeg(); + this->FillPandoraTrackIDScoreVar(propertiesMap, pfp_secondaryPCARatio, "LArThreeDPCAFeatureTool_SecondaryPCARatio"); + this->FillPandoraTrackIDScoreVar(propertiesMap, pfp_tertiaryPCARatio, "LArThreeDPCAFeatureTool_TertiaryPCARatio"); + this->FillPandoraTrackIDScoreVar(propertiesMap, pfp_vertexDist, "LArThreeDVertexDistanceFeatureTool_VertexDistance"); + } + + void Razzled::FillPandoraTrackIDScoreVar(const std::map &propertiesMap, float &var, std::string name) + { + auto propertiesMapIter = propertiesMap.find(name); + if(propertiesMapIter == propertiesMap.end()) + {}// std::cout << "Razzled Module -- Error finding variable -- " << name << std::endl; + else + var = propertiesMapIter->second; + } + + + void Razzled::FillTrackMetrics(const art::Ptr &track) + { + trk_length = track->Length(); + + const TVector3 start = track->Start(); + const TVector3 end = track->End(); + trackContained = this->InFV(start) && this->InFV(end); + + if(!fMakeTree) + return; + + trackStartX = start.X(); + trackStartY = start.Y(); + trackStartZ = start.Z(); + + trackEndX = end.X(); + trackEndY = end.Y(); + trackEndZ = end.Z(); + } + + void Razzled::FillMCSMetrics(const art::Ptr &mcs) + { + if(mcs->scatterAngles().empty()) + return; + + unsigned int counter = 0; + float maxScatter = 0.f, meanScatter = 0.f; + + for(auto const& angle : mcs->scatterAngles()) + { + if(angle < 0) + continue; + + maxScatter = std::max(maxScatter, angle); + meanScatter += angle; + counter++; + } + + if(!counter) + return; + + trk_mcsScatterMean = meanScatter / counter; + trk_mcsScatterMaxRatio = maxScatter / meanScatter; + } + + void Razzled::FillMomDiff(const art::Ptr &rangeP, const art::Ptr &mcs) + { + const float rangeMom = rangeP->range_p; + const float mcsMom = mcs->fwdMomentum(); + + trk_momDiff = (rangeMom > 0 && mcsMom > 0) ? (mcsMom - rangeMom) / rangeMom : -10.f; + } + + void Razzled::FillChi2PIDMetrics(const art::Ptr &pid) + { + const std::vector AlgScoresVec = pid->ParticleIDAlgScores(); + + std::vector> chi2s; + + for(size_t i_algscore = 0; i_algscore < AlgScoresVec.size(); i_algscore++) + { + const anab::sParticleIDAlgScores AlgScore = AlgScoresVec.at(i_algscore); + + if(AlgScore.fAlgName == "Chi2") + { + chi2s.push_back({AlgScore.fAssumedPdg, AlgScore.fValue}); + + switch(TMath::Abs(AlgScore.fAssumedPdg)) + { + case 13: + trk_chi2PIDMuon = AlgScore.fValue; + break; + case 211: + trackChi2PIDPion = AlgScore.fValue; + break; + case 321: + trackChi2PIDKaon = AlgScore.fValue; + break; + case 2212: + trk_chi2PIDProton = AlgScore.fValue; + break; + } + } + } + + if(chi2s.size() > 0) + { + std::sort(chi2s.begin(), chi2s.end(), + [](auto const& a, auto const& b) + { return a.second < b.second; }); + + chi2PDG = chi2s[0].first; + } + + trk_chi2PIDMuonPionDiff = trk_chi2PIDMuon - trackChi2PIDPion; + } + + void Razzled::FillStoppingChi2Metrics(const art::Ptr &stoppingChi2) + { + const float pol0Chi2 = stoppingChi2->pol0Chi2; + const float expChi2 = stoppingChi2->expChi2; + + trk_stoppingdEdxChi2Ratio = (pol0Chi2 > 0.f && expChi2 > 0.f) ? pol0Chi2 / expChi2 : -5.f; + trk_chi2Pol0dEdxFit = stoppingChi2->pol0Fit; + } + + void Razzled::FillDazzleMetrics(const art::Ptr &dazzle) + { + const std::map map = dazzle->mvaScoreMap; + + dazzleMuonScore = map.at(13); + dazzlePionScore = map.at(211); + dazzleProtonScore = map.at(2212); + dazzleOtherScore = map.at(0); + + dazzlePDG = dazzle->BestPDG(); + } + + void Razzled::FillShowerMetrics(const art::Ptr &shower, const std::vector> &hits) + { + const geo::GeometryCore* geom = lar::providerFrom(); + + const float length = shower->Length(); + shw_openAngle = TMath::RadToDeg() * shower->OpenAngle(); + if(shw_openAngle < 0) + shw_openAngle = -10.f; + + const TVector3 start = shower->ShowerStart(); + const TVector3 end = start + length * shower->Direction(); + + showerContained = this->InFV(start) && this->InFV(end); + + std::array showerPlaneHits = { 0, 0, 0 }; + + for(auto const& hit : hits) + showerPlaneHits[hit->WireID().Plane]++; + + std::array showerPlanePitches = { -1.f, -1.f, -1.f }; + + for(geo::PlaneGeo const& plane : geom->Iterate()) + { + const float angleToVert = geom->WireAngleToVertical(plane.View(), plane.ID()) - 0.5 * M_PI; + const float cosgamma = std::abs(std::sin(angleToVert) * shower->Direction().Y() + std::cos(angleToVert) * shower->Direction().Z()); + + showerPlanePitches[plane.ID().Plane] = plane.WirePitch() / cosgamma; + } + + int bestPlane = -1; + + if(showerPlaneHits[2] >= showerPlaneHits[1] && showerPlaneHits[2] >= showerPlaneHits[0]) + bestPlane = 2; + else if(showerPlaneHits[0] >= showerPlaneHits[1]) + bestPlane = 0; + else + bestPlane = 1; + + shw_bestdEdx = shower->dEdx()[bestPlane]; + shw_bestdEdx = std::min(shw_bestdEdx, 20.f); + shw_bestdEdx = std::max(shw_bestdEdx, -5.f); + + const float bestEnergy = shower->Energy()[bestPlane]; + const int bestPlaneHits = showerPlaneHits[bestPlane]; + const float bestPitch = showerPlanePitches[bestPlane]; + const float wiresHit = bestPitch > std::numeric_limits::epsilon() ? length / bestPitch : -5.f; + + shw_sqrtEnergyDensity = (length > 0 && bestEnergy > 0) ? std::sqrt(bestEnergy) / length : -5.f; + shw_modHitDensity = wiresHit > 1.f ? bestPlaneHits / wiresHit : -5.f; + shw_modHitDensity = std::min(shw_modHitDensity, 40.f); + + if(!fMakeTree) + return; + + showerStartX = start.X(); + showerStartY = start.Y(); + showerStartZ = start.Z(); + + showerEndX = end.X(); + showerEndY = end.Y(); + showerEndZ = end.Z(); + } + + void Razzled::FillShowerConversionGap(const art::Ptr &pfp, const std::map> &pfpMap, + const art::Ptr &shower, const art::FindOneP &pfpsToVertices) + { + const int parentId = pfp->Parent(); + auto const& parentIter = pfpMap.find(parentId); + + if(parentIter == pfpMap.end()) + return; + + const art::Ptr parent = parentIter->second; + + const art::Ptr vertex = pfpsToVertices.at(parent.key()); + + if(vertex.isNull()) + return; + + const geo::Point_t vertexPoint = vertex->position(); + const TVector3 vertexTV3 = { vertexPoint.X(), vertexPoint.Y(), vertexPoint.Z() }; + + shw_convGap = (shower->ShowerStart() - vertexTV3).Mag(); + shw_convGap = std::min(shw_convGap, 50.f); + } + + void Razzled::FillRazzleMetrics(const art::Ptr &razzle) + { + const std::map map = razzle->mvaScoreMap; + + razzleElectronScore = map.at(11); + razzlePhotonScore = map.at(22); + razzleOtherScore = map.at(0); + + razzlePDG = razzle->BestPDG(); + } + + bool Razzled::InFV(const TVector3 &pos) + { + return (pos.X() > fXMin && pos.X() < fXMax && pos.Y() > fYMin && pos.Y() < fYMax && pos.Z() > fZMin && pos.Z() < fZMax); + } + + std::map> Razzled::GetPFPMap(const std::vector> &pfps) + { + std::map> pfpMap; + for(auto const& pfp : pfps) + pfpMap[pfp->Self()] = pfp; + + return pfpMap; + } + + std::string Razzled::PDGString(const int pdg) + { + switch (std::abs(pdg)) + { + case 11: + return "Electron"; + case 13: + return "Muon"; + case 22: + return "Photon"; + case 211: + return "Pion"; + case 2212: + return "Proton"; + default: + return "Other"; + } + } + + MVAPID Razzled::RunMVA() + { + const std::vector mvaScores = reader->EvaluateMulticlass(fMethodName); + + MVAPID pidResults; + + pidResults.AddScore(11, mvaScores.at(0)); + pidResults.AddScore(13, mvaScores.at(1)); + pidResults.AddScore(22, mvaScores.at(2)); + pidResults.AddScore(211, mvaScores.at(3)); + pidResults.AddScore(2212, mvaScores.at(4)); + + if(!fMakeTree) + return pidResults; + + electronScore = mvaScores.at(0); + muonScore = mvaScores.at(1); + photonScore = mvaScores.at(2); + pionScore = mvaScores.at(3); + protonScore = mvaScores.at(4); + + bestScore = pidResults.BestScore(); + bestPDG = pidResults.BestPDG(); + + return pidResults; + } +} + +DEFINE_ART_MODULE(sbn::Razzled) diff --git a/sbncode/PID/sbn_pid.fcl b/sbncode/PID/sbn_pid.fcl index b99b0fd97..77a4a48ba 100644 --- a/sbncode/PID/sbn_pid.fcl +++ b/sbncode/PID/sbn_pid.fcl @@ -42,6 +42,33 @@ razzle: ZMax: 495 } +razzled: +{ + module_type: Razzled + PFPLabel: "pandora" + ClusterLabel: "pandora" + TrackLabel: "pandoraTrack" + ShowerLabel: "pandoraShowerSBN" + CaloLabel: "pandoraCalo" + MCSLabel: "pandoraTrackMCS" + Chi2Label: "pandoraPid" + RangeLabel: "pandoraTrackRange" + ClosestApproachLabel: "pandoraTrackClosestApproach" + StoppingChi2Label: "pandoraTrackStoppingChi2" + DazzleLabel: "pandoraTrackDazzle" + RazzleLabel: "pandoraShowerRazzle" + MinTrackLength: 0 + MinShowerEnergy: 0 + MakeTree: false + RunMVA: false + SaveFullCalo: false + XMin: -195 + XMax: 195 + YMin: -195 + YMax: 195 + ZMin: 5 + ZMax: 495 +} dazzle_sbnd: @local::dazzle dazzle_sbnd.RunMVA: true @@ -53,4 +80,17 @@ razzle_sbnd.RunMVA: true razzle_sbnd.MethodName: "BDT::BDTG" razzle_sbnd.WeightFile: "PID/Razzle.weights.xml" +razzled_sbnd: @local::razzled +razzled_sbnd.RunMVA: true +razzled_sbnd.MethodName: "BDT::BDTG" +razzled_sbnd.WeightFile: "PID/Razzled.weights.xml" + +razzled_sbnd_sce: @local::razzled_sbnd +razzled_sbnd_sce.PFPLabel: "pandoraSCE" +razzled_sbnd_sce.ClusterLabel: "pandoraSCE" +razzled_sbnd_sce.TrackLabel: "pandoraSCETrack" +razzled_sbnd_sce.ShowerLabel: "pandoraSCEShowerSBN" +razzled_sbnd_sce.CaloLabel: "pandoraSCECalo" +razzled_sbnd_sce.Chi2Label: "pandoraSCEPid" + END_PROLOG diff --git a/sbncode/PID/scripts/PFPPIDMVA.C b/sbncode/PID/scripts/PFPPIDMVA.C new file mode 100644 index 000000000..438c84726 --- /dev/null +++ b/sbncode/PID/scripts/PFPPIDMVA.C @@ -0,0 +1,88 @@ +#include "TChain.h" +#include "TFile.h" +#include "TObjString.h" +#include "TROOT.h" +#include "TString.h" +#include "TSystem.h" +#include "TTree.h" + +#include "TMVA/DataLoader.h" +#include "TMVA/Factory.h" +#include "TMVA/TMVAGui.h" +#include "TMVA/Tools.h" + +void PFPPIDMVA() +{ + gStyle->SetOptStat(0); + + TString instance = ""; + + TFile *outputFile = TFile::Open("Razzled" + instance + ".root", "RECREATE"); + + TMVA::Factory *factory = new TMVA::Factory("Razzled" + instance, outputFile, "!V:!Silent:Color:DrawProgressBar:AnalysisType=multiclass"); + TMVA::DataLoader *dataloader = new TMVA::DataLoader("Razzled" + instance); + + TChain *pfpTree = new TChain("razzled/pfpTree"); + pfpTree->Add("/ADD/WHATEVER/SAMPLES/YOU'RE/USING/HERE"); + // Recommend using a combination of a rockbox sample & an intrinsic electron neutrino sample. + + TCut generic_cuts = "!unambiguousSlice && (trk_length > 3 || showerEnergy > 10)"; + + TTree *electronTree = pfpTree->CopyTree("std::abs(truePDG)==11" + generic_cuts); + TTree *muonTree = pfpTree->CopyTree("std::abs(truePDG)==13" + generic_cuts); + TTree *photonTree = pfpTree->CopyTree("std::abs(truePDG)==22" + generic_cuts); + TTree *pionTree = pfpTree->CopyTree("std::abs(truePDG)==211" + generic_cuts); + TTree *protonTree = pfpTree->CopyTree("std::abs(truePDG)==2212" + generic_cuts); + + dataloader->AddTree(electronTree, "Electron"); + dataloader->AddTree(muonTree, "Muon"); + dataloader->AddTree(photonTree, "Photon"); + dataloader->AddTree(pionTree, "Pion"); + dataloader->AddTree(protonTree, "Proton"); + + dataloader->AddVariable("pfp_numDaughters", "PFP N Daughters", "", 'F'); + dataloader->AddVariable("pfp_maxDaughterHits", "PFP Max Daughter Hits", "", 'F'); + dataloader->AddVariable("pfp_trackScore", "PFP Track Score", "", 'F'); + + dataloader->AddVariable("trk_length", "Track Length", "cm", 'F'); + dataloader->AddVariable("trk_chi2PIDMuon", "Track Chi2 PID Muon", "", 'F'); + dataloader->AddVariable("trk_chi2PIDProton", "Track Chi2 PID Proton", "", 'F'); + dataloader->AddVariable("trk_chi2PIDMuonPionDiff", "Track Chi2 PID Muon-Pion", "", 'F'); + dataloader->AddVariable("trk_mcsScatterMean", "Track Mean MCS Scattering Angle", "mRad", 'F'); + dataloader->AddVariable("trk_mcsScatterMaxRatio", "Track Max/Mean MCS Scattering Angle Ratio", "", 'F'); + dataloader->AddVariable("trk_meanDCA", "Track Mean DCA", "cm", 'F'); + dataloader->AddVariable("trk_stoppingdEdxChi2Ratio", "Track Stopping Chi2Ratio", "", 'F'); + dataloader->AddVariable("trk_chi2Pol0dEdxFit", "Track Fitted Pol0 dE/dx", "MeV/cm", 'F'); + dataloader->AddVariable("trk_momDiff", "Track Momentum Agreement", "", 'F'); + + dataloader->AddVariable("shw_bestdEdx", "Shower Best Plane dEdx", "MeV/cm", 'F'); + dataloader->AddVariable("shw_convGap", "Shower Conversion Gap", "cm", 'F'); + dataloader->AddVariable("shw_openAngle", "Shower Opening Angle", "rad", 'F'); + dataloader->AddVariable("shw_modHitDensity", "Shower Modified Hit Density", "", 'F'); + dataloader->AddVariable("shw_sqrtEnergyDensity > 2.5 ? 2.5 : shw_sqrtEnergyDensity", "Shower Sqrt Energy Density", "", 'F'); + + const TCut baseCut("(abs(trackStartX) < 180 && abs(trackStartY) < 180 && trackStartZ > 10" + " && trackStartZ < 450 && abs(showerStartX) < 180 && abs(showerStartY) < 180" + " && showerStartZ > 10 && showerStartZ < 450 && recoPrimary == 1" + " && energyPurity > 0.5 && energyComp > 0.5)"); + + dataloader->PrepareTrainingAndTestTree(baseCut, ""); + + factory->BookMethod(dataloader, TMVA::Types::kBDT, "BDTG", + "!H:!V:NTrees=100:BoostType=Grad:Shrinkage=0.50::BaggedSampleFraction=0.60" + ":nCuts=100:MaxDepth=3:DoBoostMonitor"); + + factory->TrainAllMethods(); + + // Evaluate all MVAs using the set of test events + factory->TestAllMethods(); + + // Evaluate and compare performance of all configured MVAs + factory->EvaluateAllMethods(); + + outputFile->Close(); + + // Launch the GUI for the root macros + if (!gROOT->IsBatch()) + TMVA::TMVAMultiClassGui("Razzled" + instance + ".root"); +} diff --git a/sbncode/PID/scripts/ShowerPIDMVA.C b/sbncode/PID/scripts/ShowerPIDMVA.C index 32902c501..2b91f808b 100644 --- a/sbncode/PID/scripts/ShowerPIDMVA.C +++ b/sbncode/PID/scripts/ShowerPIDMVA.C @@ -11,54 +11,50 @@ #include "TMVA/TMVAGui.h" #include "TMVA/Tools.h" -void ShowerPIDMVA(TString fileName) +void ShowerPIDMVA() { - gStyle->SetOptStat(0); - TFile f(fileName); + gStyle->SetOptStat(0); - TFile *outputFile = TFile::Open("ShowerPIDMVA.root", "RECREATE"); + TFile *outputFile = TFile::Open("ShowerPIDMVA.root", "RECREATE"); - TMVA::Factory *factory = new TMVA::Factory("ShowerPIDMVA", outputFile, "!V:!Silent:Color:DrawProgressBar:AnalysisType=multiclass"); - TMVA::DataLoader *dataloader = new TMVA::DataLoader("dataset"); + TMVA::Factory *factory = new TMVA::Factory("ShowerPIDMVA", outputFile, "!V:!Silent:Color:DrawProgressBar:AnalysisType=multiclass"); + TMVA::DataLoader *dataloader = new TMVA::DataLoader("dataset"); - TDirectory *dir = f.GetDirectory("pandoraShowerSBNMVAPID"); - TTree *showerTree = (TTree *)dir->Get("showerTree"); + TChain *showerTree = new TChain("razzle/showerTree"); + showerTree->Add("/ADD/WHATEVER/SAMPLES/YOU'RE/USING/HERE"); - TTree *electronTree = showerTree->CopyTree("std::abs(truePdg)==11"); - TTree *photonTree = showerTree->CopyTree("std::abs(truePdg)==22"); - TTree *otherTree = showerTree->CopyTree("std::abs(truePdg)!=11 && std::abs(truePdg)!=22"); + TTree *electronTree = showerTree->CopyTree("std::abs(truePdg)==11"); + TTree *photonTree = showerTree->CopyTree("std::abs(truePdg)==22"); + TTree *otherTree = showerTree->CopyTree("std::abs(truePdg)!=11 && std::abs(truePdg)!=22"); - dataloader->AddTree(electronTree, "Electron"); - dataloader->AddTree(photonTree, "Photon"); - dataloader->AddTree(otherTree, "Other"); + dataloader->AddTree(electronTree, "Electron"); + dataloader->AddTree(photonTree, "Photon"); + dataloader->AddTree(otherTree, "Other"); - dataloader->AddVariable("bestdEdx", "Best Plane dEdx", "MeV/cm", 'F', 0, 10); - dataloader->AddVariable("convGap", "Conversion Gap", "cm", 'F', 0, 10); - dataloader->AddVariable("openAngle", "Opening Angle", "rad", 'F', 0, 10); - dataloader->AddVariable("modHitDensity", "Modified Hit Density", "", 'F', 0, 10); - dataloader->AddVariable("sqrtEnergyDensity", "Sqrt Energy Density", "", 'F', 0, 10); - // dataloader->AddVariable("logEnergyDensity", "Sqrt Energy Density", "cm^-1", 'F', 0, 10); + dataloader->AddVariable("bestdEdx", "Best Plane dEdx", "MeV/cm", 'F', 0, 10); + dataloader->AddVariable("convGap", "Conversion Gap", "cm", 'F', 0, 10); + dataloader->AddVariable("openAngle", "Opening Angle", "rad", 'F', 0, 10); + dataloader->AddVariable("modHitDensity", "Modified Hit Density", "", 'F', 0, 10); + dataloader->AddVariable("sqrtEnergyDensity", "Sqrt Energy Density", "", 'F', 0, 10); - // dataloader->AddSpectator("trueP", "True Momentum", "GeV"); - // dataloader->AddSpectator("bestEnergy", "Reco Energy", "MeV"); + const TCut baseCut("(abs(startX) < 175 && abs(startY) < 175 && startZ > 25 && startZ < 450" + "&& recoPrimary == 1 && bestEnergy>100 && trackScore < 0.5 && " + "energyPurity > 0.5 && energyComp > 0.5 && recoContained)"); - const TCut baseCut("(abs(startX) < 175 && abs(startY) < 175 && startZ > 25 && startZ < 450 && recoPrimary == 1 && bestEnergy>100 && " - "energyPurity > 0.5 && energyComp > 0.5 && recoContained)"); + dataloader->PrepareTrainingAndTestTree(baseCut, ""); - dataloader->PrepareTrainingAndTestTree(baseCut, ""); + factory->BookMethod(dataloader, TMVA::Types::kBDT, "BDTG", + "!H:!V:NTrees=100:BoostType=Grad:Shrinkage=0.50::BaggedSampleFraction=0.60:nCuts=100:MaxDepth=3:DoBoostMonitor"); - factory->BookMethod(dataloader, TMVA::Types::kBDT, "BDTG", - "!H:!V:NTrees=100:BoostType=Grad:Shrinkage=0.50::BaggedSampleFraction=0.60:nCuts=100:MaxDepth=3:DoBoostMonitor"); + factory->TrainAllMethods(); - factory->TrainAllMethods(); + factory->TestAllMethods(); - factory->TestAllMethods(); + factory->EvaluateAllMethods(); - factory->EvaluateAllMethods(); + outputFile->Close(); - outputFile->Close(); - - // Launch the GUI for the root macros - if (!gROOT->IsBatch()) - TMVA::TMVAMultiClassGui("ShowerPIDMVA.root"); + // Launch the GUI for the root macros + if (!gROOT->IsBatch()) + TMVA::TMVAMultiClassGui("ShowerPIDMVA.root"); } diff --git a/sbncode/PID/scripts/TrackPIDMVA.C b/sbncode/PID/scripts/TrackPIDMVA.C index 38b4d6b34..f3a38bde1 100644 --- a/sbncode/PID/scripts/TrackPIDMVA.C +++ b/sbncode/PID/scripts/TrackPIDMVA.C @@ -11,73 +11,61 @@ #include "TMVA/TMVAGui.h" #include "TMVA/Tools.h" -void TrackPIDMVA(TString fileName) +void TrackPIDMVA() { - gStyle->SetOptStat(0); - TFile f(fileName); + gStyle->SetOptStat(0); - TFile *outputFile = TFile::Open("TrackPIDMVA.root", "RECREATE"); + TFile *outputFile = TFile::Open("TrackPIDMVA.root", "RECREATE"); - TMVA::Factory *factory = new TMVA::Factory("TrackPIDMVA", outputFile, "!V:!Silent:Color:DrawProgressBar:AnalysisType=multiclass"); - TMVA::DataLoader *dataloader = new TMVA::DataLoader("dataset"); + TMVA::Factory *factory = new TMVA::Factory("TrackPIDMVA", outputFile, "!V:!Silent:Color:DrawProgressBar:AnalysisType=multiclass"); + TMVA::DataLoader *dataloader = new TMVA::DataLoader("dataset"); - TDirectory *dir = f.GetDirectory("pandoraTrackMVAPID"); - TTree *trackTree = (TTree *)dir->Get("trackTree"); + TChain *trackTree = new TChain("dazzle/trackTree"); + trackTree->Add("/ADD/WHATEVER/SAMPLES/YOU'RE/USING/HERE"); - TTree *muonTree = trackTree->CopyTree("std::abs(truePdg)==13"); - TTree *pionTree = trackTree->CopyTree("std::abs(truePdg)==211"); - TTree *protonTree = trackTree->CopyTree("std::abs(truePdg)==2212"); - TTree *otherTree = trackTree->CopyTree("std::abs(truePdg)!=13 && std::abs(truePdg)!=211 && std::abs(truePdg)!=2212"); + TTree *muonTree = trackTree->CopyTree("std::abs(truePdg)==13"); + TTree *pionTree = trackTree->CopyTree("std::abs(truePdg)==211"); + TTree *protonTree = trackTree->CopyTree("std::abs(truePdg)==2212"); + TTree *otherTree = trackTree->CopyTree("std::abs(truePdg)!=13 && std::abs(truePdg)!=211 && std::abs(truePdg)!=2212"); - dataloader->AddTree(muonTree, "Muon"); - dataloader->AddTree(pionTree, "Pion"); - dataloader->AddTree(protonTree, "Proton"); - dataloader->AddTree(otherTree, "Other"); + dataloader->AddTree(muonTree, "Muon"); + dataloader->AddTree(pionTree, "Pion"); + dataloader->AddTree(protonTree, "Proton"); + dataloader->AddTree(otherTree, "Other"); - dataloader->AddVariable("recoLen", "Reco. Length", "cm", 'F', 0, 250); - dataloader->AddVariable("chi2PIDMuon", "Chi2 PID Muon", "", 'F', 0, 80); - dataloader->AddVariable("chi2PIDProton", "Chi2 PID Proton", "", 'F', 0, 300); - dataloader->AddVariable("chi2PIDMuonPionDiff", "Chi2 PID Muon-Pion", "", 'F', 0, 80); + dataloader->AddVariable("recoLen", "Reco. Length", "cm", 'F', 0, 250); + dataloader->AddVariable("chi2PIDMuon", "Chi2 PID Muon", "", 'F', 0, 80); + dataloader->AddVariable("chi2PIDProton", "Chi2 PID Proton", "", 'F', 0, 300); + dataloader->AddVariable("chi2PIDMuonPionDiff", "Chi2 PID Muon-Pion", "", 'F', 0, 80); + dataloader->AddVariable("mcsScatterMean", "Mean MCS Scattering Angle", "mRad", 'F', 0, 600); + dataloader->AddVariable("mcsScatterMaxRatio", "Max/Mean MCS Scattering Angle Ratio", "", 'F', 0, 800); + dataloader->AddVariable("meanDCA", "Mean DCA", "cm", 'F', 0, 10); + dataloader->AddVariable("stoppingChi2Ratio", "Stopping CHi2Ratio", "", 'F', 0, 5); + dataloader->AddVariable("chi2Pol0Fit", "Fitted Pol0 dE/dx", "MeV/cm", 'F', 0, 20); + dataloader->AddVariable("pDiff", "Momentum Agreement", "", 'F', 0, 10); + dataloader->AddVariable("numDaughters", "Num Daughters", "", 'I', 0, 5); + dataloader->AddVariable("maxDaughterHits", "Daughter Hits", "", 'I', 0, 300); - dataloader->AddVariable("mcsScatterMean", "Mean MCS Scattering Angle", "mRad", 'F', 0, 600); - dataloader->AddVariable("mcsScatterMaxRatio", "Max/Mean MCS Scattering Angle Ratio", "", 'F', 0, 800); - dataloader->AddVariable("meanDCA", "Mean DCA", "cm", 'F', 0, 10); + const TCut baseCut("(abs(startX) < 175 && abs(startY) < 175 && startZ > 25 && startZ < 450" + "&& recoPrimary == 1 && recoLen>10 && trackScore > 0.5 && " + "energyPurity > 0.5 && energyComp > 0.5 && recoContained)"); - dataloader->AddVariable("stoppingChi2Ratio", "Stopping CHi2Ratio", "", 'F', 0, 5); - // dataloader->AddVariable("landauGaussChi2", "Landau-Gauss Chi2", "MeV/cm", 'F', 0, 300); - // dataloader->AddVariable("lgcMPV", "Landau-Gauss MPV", "MeV/cm", 'F', 0, 20); - dataloader->AddVariable("chi2Pol0Fit", "Fitted Pol0 dE/dx", "MeV/cm", 'F', 0, 20); + dataloader->PrepareTrainingAndTestTree(baseCut, ""); - dataloader->AddVariable("pDiff", "Momentum Agreement", "", 'F', 0, 10); - dataloader->AddVariable("numDaughters", "Num Daughters", "", 'I', 0, 5); - dataloader->AddVariable("maxDaughterHits", "Daughter Hits", "", 'I', 0, 300); - // dataloader->AddVariable("daughterTrackScore", "Daughter Track Score", "", 'F', 0, 300); - // dataloader->AddSpectator("trueP", "True Momentum", "GeV"); - // dataloader->AddSpectator("trueStopping", "True Stopping", ""); + factory->BookMethod(dataloader, TMVA::Types::kBDT, "BDTG", + "!H:!V:NTrees=100:BoostType=Grad:Shrinkage=0.50::BaggedSampleFraction=0.60:nCuts=100:MaxDepth=3:DoBoostMonitor"); + factory->TrainAllMethods(); - // dataloader->AddSpectator("chi2PIDPDG", "Chi2 PID Type", ""); - // dataloader->AddSpectator("chi2PIDPDGNoKaon", "Chi2 PID Type (No Kaon)", ""); + // Evaluate all MVAs using the set of test events + factory->TestAllMethods(); - const TCut baseCut("(abs(startX) < 175 && abs(startY) < 175 && startZ > 25 && startZ < 450 && recoPrimary == 1 && recoLen>0 && " - "energyPurity > 0.5 && energyComp > 0.5 && recoContained)"); + // Evaluate and compare performance of all configured MVAs + factory->EvaluateAllMethods(); - dataloader->PrepareTrainingAndTestTree(baseCut, ""); + outputFile->Close(); - - factory->BookMethod(dataloader, TMVA::Types::kBDT, "BDTG", - "!H:!V:NTrees=100:BoostType=Grad:Shrinkage=0.50::BaggedSampleFraction=0.60:nCuts=100:MaxDepth=3:DoBoostMonitor"); - factory->TrainAllMethods(); - - // Evaluate all MVAs using the set of test events - factory->TestAllMethods(); - - // Evaluate and compare performance of all configured MVAs - factory->EvaluateAllMethods(); - - outputFile->Close(); - - // Launch the GUI for the root macros - if (!gROOT->IsBatch()) - TMVA::TMVAMultiClassGui("TrackPIDMVA.root"); + // Launch the GUI for the root macros + if (!gROOT->IsBatch()) + TMVA::TMVAMultiClassGui("TrackPIDMVA.root"); }