sbn::TrackCaloSkimmer Class Reference

#include <TrackCaloSkimmer.h>

Inheritance diagram for sbn::TrackCaloSkimmer:

Classes
struct	GlobalTrackInfo

struct	Snippet

Public Member Functions
	TrackCaloSkimmer (fhicl::ParameterSet const &p)

	TrackCaloSkimmer (TrackCaloSkimmer const &)=delete

	TrackCaloSkimmer (TrackCaloSkimmer &&)=delete

TrackCaloSkimmer &	operator= (TrackCaloSkimmer const &)=delete

TrackCaloSkimmer &	operator= (TrackCaloSkimmer &&)=delete

	~TrackCaloSkimmer ()

void	analyze (art::Event const &e) override

void	respondToOpenInputFile (const art::FileBlock &fb) override

Private Member Functions
void	FillTrack (const recob::Track &track, const recob::PFParticle &pfp, float t0, const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< const recob::TrackHitMeta * > &thms, const std::vector< art::Ptr< recob::SpacePoint >> &sps, const std::vector< art::Ptr< anab::Calorimetry >> &calo, const std::map< geo::WireID, art::Ptr< raw::RawDigit >> &rawdigits, const std::vector< GlobalTrackInfo > &tracks, const geo::GeometryCore geo, const detinfo::DetectorClocksData &clock_data, const cheat::BackTrackerService bt_serv, const sbn::EDet det)

void	FillTrackDaughterRays (const recob::Track &trk, const recob::PFParticle &pfp, const std::vector< art::Ptr< recob::PFParticle >> &PFParticleList, const art::FindManyP< recob::SpacePoint > &PFParticleSPs)

void	FillTrackEndHits (const geo::GeometryCore *geometry, const detinfo::DetectorPropertiesData &dprop, const recob::Track &track, const std::vector< art::Ptr< recob::Hit >> &allHits, const art::FindManyP< recob::SpacePoint > &allHitSPs)

void	FillTrackTruth (const detinfo::DetectorClocksData &clock_data, const std::vector< art::Ptr< recob::Hit >> &trkHits, const std::vector< art::Ptr< simb::MCParticle >> &mcparticles, const std::vector< geo::BoxBoundedGeo > &active_volumes, const std::vector< std::vector< geo::BoxBoundedGeo >> &tpc_volumes, const std::map< int, std::vector< std::pair< geo::WireID, const sim::IDE * >>> id_to_ide_map, const std::map< int, std::vector< art::Ptr< recob::Hit >>> id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo)

TrackHitInfo	MakeHit (const recob::Hit &hit, unsigned hkey, const recob::TrackHitMeta &thm, const recob::Track &trk, const art::Ptr< recob::SpacePoint > &sp, const std::vector< art::Ptr< anab::Calorimetry >> &calo, const geo::GeometryCore geo, const detinfo::DetectorClocksData &dclock, const cheat::BackTrackerService bt_serv)

void	DoTailFit ()

Private Attributes
art::InputTag	fPFPproducer

std::vector< art::InputTag >	fT0producers

art::InputTag	fCALOproducer

art::InputTag	fTRKproducer

art::InputTag	fTRKHMproducer

art::InputTag	fHITproducer

std::vector< art::InputTag >	fRawDigitproducers

std::string	fG4producer

std::string	fSimChannelproducer

bool	fRequireT0

bool	fDoTailFit

bool	fVerbose

bool	fSilenceMissingDataProducts

double	fHitRawDigitsTickCollectWidth

double	fTailFitResidualRange

int	fHitRawDigitsWireCollectWidth

bool	fFillTrackEndHits

float	fTrackEndHitWireBox

float	fTrackEndHitTimeBox

std::vector< std::unique_ptr < sbn::ITCSSelectionTool > >	fSelectionTools

MetaInfo	fMeta

std::map< Snippet, int >	fSnippetCount

std::map< geo::WireID, std::pair< int, int > >	fWiresToSave

TTree *	fTree

TrackInfo *	fTrack

TFitter	fFitExp

TFitter	fFitConst

Detailed Description

Definition at line 73 of file TrackCaloSkimmer.h.

Constructor & Destructor Documentation

sbn::TrackCaloSkimmer::TrackCaloSkimmer ( fhicl::ParameterSet const & p )

explicit

Definition at line 59 of file TrackCaloSkimmer_module.cc.

   : EDAnalyzer{p},
     fFitExp(2),
     fFitConst(1)
 {
   // Grab config
   fPFPproducer  = p.get< art::InputTag > ("PFPproducer","pandoraGausCryo0");
   fT0producers   = p.get< std::vector<art::InputTag> > ("T0producers", {"pandoraGausCryo0"} );
   fCALOproducer = p.get< art::InputTag > ("CALOproducer");
   fTRKproducer  = p.get< art::InputTag > ("TRKproducer" );
   fTRKHMproducer= p.get< art::InputTag   > ("TRKHMproducer", "");
   fHITproducer  = p.get< art::InputTag > ("HITproducer" );
   fG4producer  = p.get< std::string > ("G4producer" );
   fSimChannelproducer  = p.get< std::string > ("SimChannelproducer" );
   fRequireT0 = p.get<bool>("RequireT0", false);
   fDoTailFit = p.get<bool>("DoTailFit", true);
   fVerbose = p.get<bool>("Verbose", false);
   fSilenceMissingDataProducts = p.get<bool>("SilenceMissingDataProducts", false);
   fHitRawDigitsTickCollectWidth = p.get<double>("HitRawDigitsTickCollectWidth", 50.);
   fHitRawDigitsWireCollectWidth = p.get<int>("HitRawDigitsWireCollectWidth", 5);
   fTailFitResidualRange = p.get<double>("TailFitResidualRange", 5.);
   if (fTailFitResidualRange > 10.) {
     std::cout << "sbn::TrackCaloSkimmer: Bad tail fit residual range config :(" << fTailFitResidualRange << "). Fits will not be meaningful.\n";
   }
   fFillTrackEndHits = p.get<bool>("FillTrackEndHits", true);
   fTrackEndHitWireBox = p.get<float>("TrackEndHitWireBox", 60); // 20 cm
   fTrackEndHitTimeBox = p.get<float>("TrackEndHitTimeBox", 300); // about 20cm
 
   fRawDigitproducers = p.get<std::vector<art::InputTag>>("RawDigitproducers", {});
 
   std::vector<fhicl::ParameterSet> selection_tool_configs(p.get<std::vector<fhicl::ParameterSet>>("SelectionTools"), {});
   for (const fhicl::ParameterSet &p: selection_tool_configs) {
     fSelectionTools.push_back(art::make_tool<sbn::ITCSSelectionTool>(p));
   }
 
   // Setup meta info
   fMeta.iproc = -1;
   fMeta.ifile = -1;
   const char *process_str = std::getenv("PROCESS");
   if (process_str) {
     try {
       fMeta.iproc = std::stoi(process_str);
     }
     catch (...) {}
   }
 
   // Output stuff
   fTrack = new sbn::TrackInfo();
 
   art::ServiceHandle<art::TFileService> tfs;
   fTree = tfs->make<TTree>("TrackCaloSkim", "Calo Tree");
   fTree->Branch("trk", &fTrack);
 }

sbn::TrackCaloSkimmer::TrackCaloSkimmer ( TrackCaloSkimmer const & )

delete

sbn::TrackCaloSkimmer::TrackCaloSkimmer ( TrackCaloSkimmer && )

delete

sbn::TrackCaloSkimmer::~TrackCaloSkimmer ( )

Definition at line 55 of file TrackCaloSkimmer_module.cc.

                                        {
   delete fTrack;
 }

Member Function Documentation

void sbn::TrackCaloSkimmer::analyze ( art::Event const & e )

override

Definition at line 113 of file TrackCaloSkimmer_module.cc.

 {
   unsigned evt = e.event();
   unsigned sub = e.subRun();
   unsigned run = e.run();
   if (fVerbose) {
     std::cout << "[TrackCaloSkimmer::analyzeEvent] Run: " << run << ", SubRun: " << sub << ", Event: "<< evt << ", Is Data: " << e.isRealData() << std::endl;
   }
 
   fMeta.evt = evt;
   fMeta.subrun = sub;
   fMeta.run = run;
   fMeta.time = e.time().value();
 
   // Services
   const geo::GeometryCore *geometry = lar::providerFrom<geo::Geometry>();
   auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
   auto const dprop =
     art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e, clock_data);
 
   // Identify which detector: can only detect either sbnd or icarus
 
   std::string gdml = geometry->GDMLFile();
   gdml = basename(gdml.c_str()); 
 
   for(unsigned int i = 0; i <gdml.size(); ++i) gdml[i] = std::tolower(gdml[i]); 
 
   EDet det = kNOTDEFINED;
 
   const bool hasSBND = ((gdml.find("sbnd") != std::string::npos) ||
                         (geometry->DetectorName().find("sbnd") != std::string::npos));
 
   const bool hasICARUS = ((gdml.find("icarus") != std::string::npos) ||
                         (geometry->DetectorName().find("icarus") != std::string::npos));
 
   if(hasSBND == hasICARUS) { 
     std::cout << "TrackCaloSkimmer: Unable to automatically determine either SBND or ICARUS!" << std::endl;
     abort();
   }
  
   if(hasSBND) det = kSBND;
   if(hasICARUS) det = kICARUS;
 
   // Setup the volumes
   std::vector<std::vector<geo::BoxBoundedGeo>> TPCVols;
   std::vector<geo::BoxBoundedGeo> AVs;
 
   // First the TPC
   for (auto const &cryo: geometry->IterateCryostats()) {
     geo::GeometryCore::TPC_iterator iTPC = geometry->begin_TPC(cryo.ID()),
                                     tend = geometry->end_TPC(cryo.ID());
     std::vector<geo::BoxBoundedGeo> this_tpc_volumes;
     while (iTPC != tend) {
       geo::TPCGeo const& TPC = *iTPC;
       this_tpc_volumes.push_back(TPC.ActiveBoundingBox());
       iTPC++;
     }
      TPCVols.push_back(std::move(this_tpc_volumes));
   }
 
   for (const std::vector<geo::BoxBoundedGeo> &tpcs: TPCVols) {
     double XMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinX() < rhs.MinX(); })->MinX();
     double YMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinY() < rhs.MinY(); })->MinY();
     double ZMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinZ() < rhs.MinZ(); })->MinZ();
 
     double XMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxX() < rhs.MaxX(); })->MaxX();
     double YMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxY() < rhs.MaxY(); })->MaxY();
     double ZMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxZ() < rhs.MaxZ(); })->MaxZ();
 
     AVs.emplace_back(XMin, XMax, YMin, YMax, ZMin, ZMax);
   }
 
   // Truth information
   std::vector<art::Ptr<simb::MCParticle>> mcparticles;
   if (fG4producer.size()) {
     art::ValidHandle<std::vector<simb::MCParticle>> mcparticle_handle = e.getValidHandle<std::vector<simb::MCParticle>>(fG4producer);
     art::fill_ptr_vector(mcparticles, mcparticle_handle);
   }
 
   std::vector<art::Ptr<sim::SimChannel>> simchannels;
   if (fSimChannelproducer.size()) {
     art::ValidHandle<std::vector<sim::SimChannel>> simchannel_handle = e.getValidHandle<std::vector<sim::SimChannel>>(fSimChannelproducer);
     art::fill_ptr_vector(simchannels, simchannel_handle);
   }
 
   // Reconstructed Information
   std::vector<art::Ptr<recob::PFParticle>> PFParticleList;
   try {
     art::ValidHandle<std::vector<recob::PFParticle>> pfparticles = e.getValidHandle<std::vector<recob::PFParticle>>(fPFPproducer);
     art::fill_ptr_vector(PFParticleList, pfparticles);
   }
   catch(...) {
       std::cout << "PFP's with tag: " << fPFPproducer << " not present.\n";
       // Real data may have missing products -- just ignore the event
       if (fSilenceMissingDataProducts) return;
       else throw;
   }
 
   // PFP-associated data
   std::vector<art::FindManyP<anab::T0>> fmT0;
   for (unsigned i_label = 0; i_label < fT0producers.size(); i_label++) {
     fmT0.emplace_back(PFParticleList, e, fT0producers[i_label]);
   }
   art::FindManyP<recob::SpacePoint> PFParticleSPs(PFParticleList, e, fPFPproducer);
 
   // Now we don't need to guard access to further data. If this is an empty event it should be caught by PFP's or Hit's
   art::ValidHandle<std::vector<recob::Track>> tracks = e.getValidHandle<std::vector<recob::Track>>(fTRKproducer); 
 
   // Track - associated data
   art::FindManyP<recob::Track> fmTracks(PFParticleList, e, fTRKproducer);
 
   art::InputTag thm_label = fTRKHMproducer.empty() ? fTRKproducer : fTRKHMproducer;
   art::FindManyP<recob::Hit, recob::TrackHitMeta> fmtrkHits(tracks, e, thm_label);
   art::FindManyP<anab::Calorimetry> fmCalo(tracks, e, fCALOproducer);
 
   // Collect raw digits for saving hits
   std::vector<art::Ptr<raw::RawDigit>> rawdigitlist;
   for (const art::InputTag &t: fRawDigitproducers) {
     art::ValidHandle<std::vector<raw::RawDigit>> thisdigits = e.getValidHandle<std::vector<raw::RawDigit>>(t);
     art::fill_ptr_vector(rawdigitlist, thisdigits);
   }
 
   // The raw digit list is not sorted, so make it into a map on the WireID
   std::map<geo::WireID, art::Ptr<raw::RawDigit>> rawdigits;
   for (const art::Ptr<raw::RawDigit> &d: rawdigitlist) {
     
     std::vector<geo::WireID> wids;
     // Handle bad channel ID
     try {
       wids = geometry->ChannelToWire(d->Channel());
     }
     catch(...) {
       continue;
     }
 
     // Ignore channel with no mapped wire
     if (wids.size() == 0) continue;
 
     // Ignore wires that are already mapped
     if (rawdigits.count(wids[0])) continue;
 
     rawdigits[wids[0]] = d;
   } 
 
   // Collect all hits
   art::ValidHandle<std::vector<recob::Hit>> allhit_handle = e.getValidHandle<std::vector<recob::Hit>>(fHITproducer);
   std::vector<art::Ptr<recob::Hit>> allHits;
   art::fill_ptr_vector(allHits, allhit_handle);
 
   // And lookup the SP's
   art::FindManyP<recob::SpacePoint> allHitSPs(allHits, e, fPFPproducer);
 
   // Prep truth-to-reco-matching info
   //
   // Use helper functions from CAFMaker/FillTrue
   std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> id_to_ide_map;
   std::map<int, std::vector<art::Ptr<recob::Hit>>> id_to_truehit_map;
   const cheat::BackTrackerService *bt = NULL;
 
   if (simchannels.size()) {
     art::ServiceHandle<cheat::BackTrackerService> bt_serv;
     id_to_ide_map = caf::PrepSimChannels(simchannels, *geometry);
     id_to_truehit_map = caf::PrepTrueHits(allHits, clock_data, *bt_serv.get());
     bt = bt_serv.get();
   }
 
   // service data
 
   // Build global track info
   std::vector<GlobalTrackInfo> track_infos;
   for (const recob::Track &t: *tracks) {
     track_infos.push_back({
       t.Start(), t.End(), t.StartDirection(), t.EndDirection(), t.ID()
     });
   }
 
   for (art::Ptr<recob::PFParticle> p_pfp: PFParticleList) {
     const recob::PFParticle &pfp = *p_pfp;
 
     const std::vector<art::Ptr<recob::Track>> thisTrack = fmTracks.at(p_pfp.key());
     if (thisTrack.size() != 1)
       continue;
 
     art::Ptr<recob::Track> trkPtr = thisTrack.at(0);
 
     std::vector<art::Ptr<anab::Calorimetry>> emptyCaloVector;
     const std::vector<art::Ptr<anab::Calorimetry>> &calo = fmCalo.isValid() ? fmCalo.at(trkPtr.key()) : emptyCaloVector;
 
     std::vector<art::Ptr<recob::Hit>> emptyHitVector;
     const std::vector<art::Ptr<recob::Hit>> &trkHits  = fmtrkHits.isValid() ? fmtrkHits.at(trkPtr.key()) : emptyHitVector;
 
     art::FindManyP<recob::SpacePoint> fmtrkHitSPs(trkHits, e, fPFPproducer);
 
     std::vector<const recob::TrackHitMeta*> emptyTHMVector;
     const std::vector<const recob::TrackHitMeta*> &trkHitMetas = fmtrkHits.isValid() ? fmtrkHits.data(trkPtr.key()) : emptyTHMVector;
 
     art::Ptr<recob::SpacePoint> nullSP;
     std::vector<art::Ptr<recob::SpacePoint>> trkHitSPs;
     if (fmtrkHitSPs.isValid()) {
       for (unsigned i_hit = 0; i_hit < trkHits.size(); i_hit++) {
         const std::vector<art::Ptr<recob::SpacePoint>> &h_sp = fmtrkHitSPs.at(i_hit);
         if (h_sp.size()) {
           trkHitSPs.push_back(h_sp.at(0));
         }
         else {
           trkHitSPs.push_back(nullSP);
         }
       }
     }
 
     int whicht0 = -1;
     float t0 = std::numeric_limits<float>::signaling_NaN();
     for (unsigned i_t0 = 0; i_t0 < fmT0.size(); i_t0++) {
       if (fmT0[i_t0].isValid() && fmT0[i_t0].at(p_pfp.key()).size()) {
         t0 = fmT0[i_t0].at(p_pfp.key()).at(0)->Time();
         whicht0 = i_t0;
 
         if (fVerbose) std::cout << "Track: " << trkPtr->ID() << " Has T0 (" << fT0producers[i_t0] << ")\n";
 
         break;
       }
     }
 
     if (fRequireT0 && whicht0 < 0) {
       continue;
     }
 
     if (fVerbose) std::cout << "Processing new track! ID: " << trkPtr->ID() << " time: " << t0 << std::endl;
 
     // Reset the track object
     *fTrack = sbn::TrackInfo();
 
     // Reset other persistent info
     fSnippetCount.clear();
     fWiresToSave.clear();
 
     // Fill the track!
     FillTrack(*trkPtr, pfp, t0, trkHits, trkHitMetas, trkHitSPs, calo, rawdigits, track_infos, geometry, clock_data, bt, det);
     fTrack->whicht0 = whicht0;
 
     FillTrackDaughterRays(*trkPtr, pfp, PFParticleList, PFParticleSPs);
 
     if (fFillTrackEndHits) FillTrackEndHits(geometry, dprop, *trkPtr, allHits, allHitSPs);
 
     // Fill the truth information if configured
     if (simchannels.size()) FillTrackTruth(clock_data, trkHits, mcparticles, AVs, TPCVols, id_to_ide_map, id_to_truehit_map, dprop, geometry); 
 
     // Save?
     bool select = false;
     if (!fSelectionTools.size()) select = true;
 
     // Take the OR of each selection tool
     int i_select = 0;
     for (const std::unique_ptr<sbn::ITCSSelectionTool> &t: fSelectionTools) {
       if (t->DoSelect(*fTrack)) {
         select = true;
         fTrack->selected = i_select;
         fTrack->nprescale = t->GetPrescale();
         break;
       }
       i_select ++;
     }
 
     // Save!
     if (select) {
       if (fVerbose) std::cout << "Track Selected! By tool: " << i_select << std::endl;
       fTree->Fill();
     }
   }
 }

void sbn::TrackCaloSkimmer::DoTailFit ( )

private

Definition at line 1114 of file TrackCaloSkimmer_module.cc.

                                     {
   // Try fitting the constant and exponentials to the tail of dQ/dx v. RR on the collection plane
   std::vector<double> fit_rr;
   std::vector<double> fit_dqdx;
 
   for (const TrackHitInfo &h: fTrack->hits2) {
     if (h.oncalo && h.rr > 0. && h.rr < fTailFitResidualRange) {
       fit_rr.push_back(h.rr);
       fit_dqdx.push_back(h.dqdx);
     }
   }
 
   // TODO: should we throw an exception here??
   // 
   // If there is too much data to fit in the array, throw exception
   if (fit_rr.size() > MAX_N_FIT_DATA) {
     throw cet::exception("sbn::TrackCaloSkimmer::DoTailFit: More fitting points required (" 
       + std::to_string(fit_rr.size()) + ") than available in fit array (" + std::to_string(MAX_N_FIT_DATA) + ").\n");
   }
 
   // Copy the fit data to the global array
   for (unsigned i = 0; i < fit_rr.size() && i < MAX_N_FIT_DATA; i++) {
     FIT_RR[i] = fit_rr[i];
     FIT_DQDX[i] = fit_dqdx[i];
   }
   N_FIT_DATA = std::min(fit_rr.size(), MAX_N_FIT_DATA);
 
   fTrack->n_fit_point = N_FIT_DATA;
   if (fTrack->n_fit_point > 2) { // need more points than params
     // Fit the Exponential
     fFitExp.SetFCN(ExpResiduals);
     fFitExp.SetParameter(0, "A", *std::max_element(fit_dqdx.begin(), fit_dqdx.end()), 200, 0, 5000);
     fFitExp.SetParameter(1, "R", 10., 0.5, 0, 1000);
     fFitExp.ExecuteCommand("MIGRAD", 0, 0);
 
     double A = fFitExp.GetParameter(0);
     double R = fFitExp.GetParameter(1);
 
     int nparam;
     double param[2] {A, R};
     double residuals = -1;
     ExpResiduals(nparam, NULL, residuals, param, 0);
 
     fTrack->exp_fit_A = A; 
     fTrack->exp_fit_R = R;
     fTrack->exp_fit_residuals = residuals;
 
     // Fit the Constant
     fFitConst.SetFCN(ConstResiduals);
     fFitConst.SetParameter(0, "C", std::accumulate(fit_dqdx.begin(), fit_dqdx.end(), 0.), 200, 0, 5000);
     fFitConst.ExecuteCommand("MIGRAD", 0, 0);
 
     double C = fFitConst.GetParameter(0);
 
     double cresiduals = -1;
     ConstResiduals(nparam, NULL, cresiduals, &C, 0);
 
     fTrack->const_fit_C = C;    
     fTrack->const_fit_residuals = cresiduals;
   }
 }

void sbn::TrackCaloSkimmer::FillTrack	(	const recob::Track &	track,
		const recob::PFParticle &	pfp,
		float	t0,
		const std::vector< art::Ptr< recob::Hit >> &	hits,
		const std::vector< const recob::TrackHitMeta * > &	thms,
		const std::vector< art::Ptr< recob::SpacePoint >> &	sps,
		const std::vector< art::Ptr< anab::Calorimetry >> &	calo,
		const std::map< geo::WireID, art::Ptr< raw::RawDigit >> &	rawdigits,
		const std::vector< GlobalTrackInfo > &	tracks,
		const geo::GeometryCore *	geo,
		const detinfo::DetectorClocksData &	clock_data,
		const cheat::BackTrackerService *	bt_serv,
		const sbn::EDet	det
	)

private

Definition at line 985 of file TrackCaloSkimmer_module.cc.

                        {
 
   // Fill top level stuff
   fTrack->meta = fMeta;
   fTrack->t0 = t0;
   fTrack->id = track.ID();
   fTrack->clear_cosmic_muon = pfp.Parent() == recob::PFParticle::kPFParticlePrimary;
 
   fTrack->length = track.Length();
   fTrack->start.x = track.Start().X();
   fTrack->start.y = track.Start().Y();
   fTrack->start.z = track.Start().Z();
   fTrack->end.x = track.End().X();
   fTrack->end.y = track.End().Y();
   fTrack->end.z = track.End().Z();
   fTrack->dir.x = track.StartDirection().X();
   fTrack->dir.y = track.StartDirection().Y();
   fTrack->dir.z = track.StartDirection().Z();
 
   if (hits.size() > 0) {
     fTrack->cryostat = hits[0]->WireID().Cryostat;
   }
 
   // Fill each hit
   for (unsigned i_hit = 0; i_hit < hits.size(); i_hit++) {
     sbn::TrackHitInfo hinfo = MakeHit(*hits[i_hit], hits[i_hit].key(), *thms[i_hit], track, sps[i_hit], calo, geo, clock_data, bt_serv);
     if (hinfo.h.plane == 0) {
       fTrack->hits0.push_back(hinfo);
     }
     else if (hinfo.h.plane == 1) {
       fTrack->hits1.push_back(hinfo);
     } 
     else if (hinfo.h.plane == 2) {
       fTrack->hits2.push_back(hinfo);
     }
   }
 
   // Hit summary info
   fTrack->hit_min_time_p0_tpcE = HitMinTime(fTrack->hits0, true, det);
   fTrack->hit_max_time_p0_tpcE = HitMaxTime(fTrack->hits0, true, det);
   fTrack->hit_min_time_p0_tpcW = HitMinTime(fTrack->hits0, false, det);
   fTrack->hit_max_time_p0_tpcW = HitMaxTime(fTrack->hits0, false, det);
   fTrack->hit_min_time_p1_tpcE = HitMinTime(fTrack->hits1, true, det);
   fTrack->hit_max_time_p1_tpcE = HitMaxTime(fTrack->hits1, true, det);
   fTrack->hit_min_time_p1_tpcW = HitMinTime(fTrack->hits1, false, det);
   fTrack->hit_max_time_p1_tpcW = HitMaxTime(fTrack->hits1, false, det);
   fTrack->hit_min_time_p2_tpcE = HitMinTime(fTrack->hits2, true, det);
   fTrack->hit_max_time_p2_tpcE = HitMaxTime(fTrack->hits2, true, det);
   fTrack->hit_min_time_p2_tpcW = HitMinTime(fTrack->hits2, false, det);
   fTrack->hit_max_time_p2_tpcW = HitMaxTime(fTrack->hits2, false, det);
 
   // Save information on a fit to the end of the track
   if (fDoTailFit) DoTailFit();
 
   // Save the Wire ADC values we need to
   for (auto const &w_pair: fWiresToSave) {
     geo::WireID wire = w_pair.first;
 
     if (rawdigits.count(wire)) {
       const raw::RawDigit &thisdigit = *rawdigits.at(wire);
       int min_tick = std::max(0, w_pair.second.first);
       int max_tick = std::min((int)thisdigit.NADC(), w_pair.second.second);
 
       // collect the adcs
       std::vector<short> adcs;
       for (int t = min_tick; t < max_tick; t++) {
         adcs.push_back(thisdigit.ADC(t));
       }
 
       WireInfo winfo;
       winfo.wire = wire.Wire;
       winfo.plane = wire.Plane;
       winfo.tpc = wire.TPC;
       winfo.channel = geo->PlaneWireToChannel(wire);
       winfo.tdc0 = min_tick;
       winfo.adcs = adcs;
 
       if (winfo.plane == 0) {
         fTrack->wires0.push_back(winfo);
       }
       else if (winfo.plane == 1) {
         fTrack->wires1.push_back(winfo);
       }
       else if (winfo.plane == 2) {
         fTrack->wires2.push_back(winfo);
       }
     }
   }
   // Sort the ADC values by wire
   std::sort(fTrack->wires0.begin(), fTrack->wires0.end(), [](auto const &lhs, auto const &rhs) {return lhs.wire < rhs.wire;});
   std::sort(fTrack->wires1.begin(), fTrack->wires1.end(), [](auto const &lhs, auto const &rhs) {return lhs.wire < rhs.wire;});
   std::sort(fTrack->wires2.begin(), fTrack->wires2.end(), [](auto const &lhs, auto const &rhs) {return lhs.wire < rhs.wire;});
 
   // get information on nearby tracks
   for (const GlobalTrackInfo &othr: tracks) {
     if (othr.ID == track.ID()) continue;
 
     if ((track.End() - othr.start).r() < 50. || (track.End() - othr.end).r() < 50.) {
       fTrack->tracks_near_end_dist.push_back(std::min((track.End() - othr.start).r(), (track.End() - othr.end).r())); 
       fTrack->tracks_near_end_costh.push_back(
         (track.End() - othr.start).r() < (track.End() - othr.end).r() ? 
           track.EndDirection().Dot(othr.dir) : track.EndDirection().Dot(othr.enddir)); 
     }
   }
 
   for (const GlobalTrackInfo &othr: tracks) {
     if (othr.ID == track.ID()) continue;
 
     if ((track.Start() - othr.start).r() < 50. || (track.Start() - othr.end).r() < 50.) {
       fTrack->tracks_near_start_dist.push_back(std::min((track.Start() - othr.start).r(), (track.Start() - othr.end).r())); 
       fTrack->tracks_near_start_costh.push_back(
         (track.Start() - othr.start).r() < (track.Start() - othr.end).r() ? 
           track.StartDirection().Dot(othr.dir) : track.StartDirection().Dot(othr.enddir)); 
     }
   }
 
 }

void sbn::TrackCaloSkimmer::FillTrackDaughterRays	(	const recob::Track &	trk,
		const recob::PFParticle &	pfp,
		const std::vector< art::Ptr< recob::PFParticle >> &	PFParticleList,
		const art::FindManyP< recob::SpacePoint > &	PFParticleSPs
	)

private

Definition at line 960 of file TrackCaloSkimmer_module.cc.

                                                         {
 
   for (unsigned d: pfp.Daughters()) {
     const recob::PFParticle &d_pfp = *PFParticleList[d];
 
     fTrack->daughter_pdg.push_back(d_pfp.PdgCode());
 
     unsigned nsp = 0;
     float min_distance = -1.;
     for (const art::Ptr<recob::SpacePoint> &sp: PFParticleSPs.at(d)) {
       if (min_distance < 0. || (sp->position() - trk.End()).r() < min_distance) {
         min_distance = (sp->position() - trk.End()).r();
       }
       nsp++;
     }
 
     fTrack->daughter_sp_toend_dist.push_back(min_distance);
     fTrack->daughter_nsp.push_back(nsp);
   }
 
 }

void sbn::TrackCaloSkimmer::FillTrackEndHits	(	const geo::GeometryCore *	geometry,
		const detinfo::DetectorPropertiesData &	dprop,
		const recob::Track &	track,
		const std::vector< art::Ptr< recob::Hit >> &	allHits,
		const art::FindManyP< recob::SpacePoint > &	allHitSPs
	)

private

Definition at line 835 of file TrackCaloSkimmer_module.cc.

                                                     {
 
   (void) dprop; // TODO: use??
 
   geo::TPCID tpc_end = geometry->FindTPCAtPosition(track.End());
   if (!tpc_end) return;
 
   geo::PlaneID plane_end(tpc_end, 2 /* collection */);
 
   float end_w = geometry->WireCoordinate(track.End(), plane_end);
 
   float end_t = -1000.;
   float closest_wire_dist = -1.;
   // Get the hit closest to the end to get the end time 
   for (const TrackHitInfo &h: fTrack->hits2) {
     if (h.oncalo && (closest_wire_dist < 0. || abs(h.h.wire - end_w) < closest_wire_dist)) {
       closest_wire_dist = abs(h.h.wire - end_w);
       end_t = h.h.time;
     }
   }
 
   for (const art::Ptr<recob::Hit> &hit: allHits) {
     geo::PlaneID h_p = hit->WireID();
     if (h_p != plane_end) continue;
 
     // Inside the box?
     float h_w = (float)hit->WireID().Wire;
     float h_t = hit->PeakTime();
 
     if (abs(h_w - end_w) < fTrackEndHitWireBox &&
         abs(h_t - end_t) < fTrackEndHitTimeBox) {
 
       // HitInfo to save
       sbn::HitInfo hinfo;
       
       // information from the hit object
       hinfo.integral = hit->Integral();
       hinfo.sumadc = hit->SummedADC();
       hinfo.width = hit->RMS();
       hinfo.time = hit->PeakTime();
       hinfo.mult = hit->Multiplicity();
       hinfo.wire = hit->WireID().Wire;
       hinfo.plane = hit->WireID().Plane;
       hinfo.tpc = hit->WireID().TPC;
       hinfo.end = hit->EndTick();
       hinfo.start = hit->StartTick();
       hinfo.id = hit.key();
 
       const std::vector<art::Ptr<recob::SpacePoint>> &h_sp = allHitSPs.at(hit.key());
       if (h_sp.size()) {
         const recob::SpacePoint &sp = *h_sp[0];
         hinfo.sp.x = sp.position().x();
         hinfo.sp.y = sp.position().y();
         hinfo.sp.z = sp.position().z();
 
         hinfo.hasSP = true;
       }
       else {
         hinfo.hasSP = false;
       }
 
       fTrack->endhits.push_back(hinfo);
       
     } 
   }
 
 }

void sbn::TrackCaloSkimmer::FillTrackTruth	(	const detinfo::DetectorClocksData &	clock_data,
		const std::vector< art::Ptr< recob::Hit >> &	trkHits,
		const std::vector< art::Ptr< simb::MCParticle >> &	mcparticles,
		const std::vector< geo::BoxBoundedGeo > &	active_volumes,
		const std::vector< std::vector< geo::BoxBoundedGeo >> &	tpc_volumes,
		const std::map< int, std::vector< std::pair< geo::WireID, const sim::IDE * >>>	id_to_ide_map,
		const std::map< int, std::vector< art::Ptr< recob::Hit >>>	id_to_truehit_map,
		const detinfo::DetectorPropertiesData &	dprop,
		const geo::GeometryCore *	geo
	)

private

Definition at line 907 of file TrackCaloSkimmer_module.cc.

                                 {
 
   // Lookup the true-particle match -- use utils in CAF
   std::vector<std::pair<int, float>> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clock_data, trkHits, true);
   float total_energy = CAFRecoUtils::TotalHitEnergy(clock_data, trkHits);
 
   fTrack->truth.depE = total_energy / 1000. /* MeV -> GeV */;
   
   // sort highest energy match to lowest
   std::sort(matches.begin(), matches.end(),
       [](const auto &a, const auto &b) {
         return a.second > b.second;
       }
   );
 
   // Save the best match
   if (matches.size()) {
     std::pair<int, float> bestmatch = matches[0];
 
     fTrack->truth.pur = bestmatch.second / total_energy;
 
     for (const art::Ptr<simb::MCParticle> &p_mcp: mcparticles) {
       if (p_mcp->TrackId() == bestmatch.first) {
         if (fVerbose) std::cout << "Matched! Track ID: " << p_mcp->TrackId() << " pdg: " << p_mcp->PdgCode() << " process: " << p_mcp->EndProcess() << std::endl;
         fTrack->truth.p = TrueParticleInfo(*p_mcp, active_volumes, tpc_volumes, id_to_ide_map, id_to_truehit_map, dprop, geo);
         fTrack->truth.eff = fTrack->truth.depE / (fTrack->truth.p.plane0VisE + fTrack->truth.p.plane1VisE + fTrack->truth.p.plane2VisE);
 
         // Lookup any Michel
         for (const art::Ptr<simb::MCParticle> &d_mcp: mcparticles) {
           if (d_mcp->Mother() == p_mcp->TrackId() && // correct parent
               (d_mcp->Process() == "Decay" || d_mcp->Process() == "muMinusCaptureAtRest") && // correct process
               abs(d_mcp->PdgCode()) == 11) { // correct PDG code
 
             fTrack->truth.michel = TrueParticleInfo(*d_mcp, active_volumes, tpc_volumes, id_to_ide_map, id_to_truehit_map, dprop, geo);
             break;
           }
         }
 
         break;
       } 
     }
   }
 }

sbn::TrackHitInfo sbn::TrackCaloSkimmer::MakeHit	(	const recob::Hit &	hit,
		unsigned	hkey,
		const recob::TrackHitMeta &	thm,
		const recob::Track &	trk,
		const art::Ptr< recob::SpacePoint > &	sp,
		const std::vector< art::Ptr< anab::Calorimetry >> &	calo,
		const geo::GeometryCore *	geo,
		const detinfo::DetectorClocksData &	dclock,
		const cheat::BackTrackerService *	bt_serv
	)

private

Definition at line 1177 of file TrackCaloSkimmer_module.cc.

                                             {
 
   // TrackHitInfo to save
   sbn::TrackHitInfo hinfo;
 
   // information from the hit object
   hinfo.h.integral = hit.Integral();
   hinfo.h.sumadc = hit.SummedADC();
   hinfo.h.width = hit.RMS();
   hinfo.h.time = hit.PeakTime();
   hinfo.h.mult = hit.Multiplicity();
   hinfo.h.wire = hit.WireID().Wire;
   hinfo.h.plane = hit.WireID().Plane;
   hinfo.h.channel = geo->PlaneWireToChannel(hit.WireID());
   hinfo.h.tpc = hit.WireID().TPC;
   hinfo.h.end = hit.EndTick();
   hinfo.h.start = hit.StartTick();
   hinfo.h.id = (int)hkey;
 
   // Do back-tracking on each hit
   if (bt_serv) {
     // The default BackTracking function goes from (peak - width, peak + width).
     //
     // This time range does not match well hits with a non-Gaussian shape where
     // the Gaussian-fit-width does not replicate the width of the pulse. 
     //
     // Instead, we use the Hit (start, end) time range. This is also more relevant
     // for (e.g.) the SummedADC charge extraction method.
     //
     // Don't use this:
     // std::vector<sim::TrackIDE> ides = bt_serv->HitToTrackIDEs(dclock, hit);
     //
     // Use this:
     std::vector<sim::TrackIDE> ides = bt_serv->ChannelToTrackIDEs(dclock, hit.Channel(), hit.StartTick(), hit.EndTick());
 
     hinfo.h.truth.e = 0.;
     hinfo.h.truth.nelec = 0.;
 
     for (const sim::TrackIDE &ide: ides) {
       hinfo.h.truth.e += ide.energy;
       hinfo.h.truth.nelec += ide.numElectrons;
     }
   }
   else {
     hinfo.h.truth.e = -1.;
     hinfo.h.truth.nelec = -1.;
   }
 
   // look up the snippet
   sbn::TrackCaloSkimmer::Snippet snippet {hit.WireID(), hit.StartTick(), hit.EndTick()};
   if (!fSnippetCount.count(snippet)) {
     fSnippetCount[snippet] = 0;
     hinfo.i_snippet = 0;
   }
   else {
     fSnippetCount[snippet] ++;
     hinfo.i_snippet = fSnippetCount[snippet];
   }
 
   // Which wires to save
   int min_tick = (int)std::floor(hit.PeakTime() - fHitRawDigitsTickCollectWidth);
   int max_tick = (int)std::ceil(hit.PeakTime() + fHitRawDigitsTickCollectWidth);
   for (int wire = hinfo.h.wire - fHitRawDigitsWireCollectWidth; wire <= hinfo.h.wire + fHitRawDigitsWireCollectWidth; wire++) {
     geo::WireID w(hit.WireID(), wire);
 
     if (fWiresToSave.count(w)) {
       fWiresToSave.at(w).first = std::min(fWiresToSave.at(w).first, min_tick);
       fWiresToSave.at(w).second = std::max(fWiresToSave.at(w).second, max_tick);
     }
     else {
       fWiresToSave[w] = {min_tick, max_tick};
     }
   }
 
   // Information from the TrackHitMeta
   bool badhit = (thm.Index() == std::numeric_limits<unsigned int>::max()) ||
                     (!trk.HasValidPoint(thm.Index()));
 
   hinfo.ontraj = !badhit;
 
   // Save trajectory information if we can
   if (!badhit) {
     geo::Point_t loc = trk.LocationAtPoint(thm.Index());
     hinfo.tp.x = loc.X();
     hinfo.tp.y = loc.Y();
     hinfo.tp.z = loc.Z();
 
     geo::Vector_t dir = trk.DirectionAtPoint(thm.Index());
     hinfo.dir.x = dir.X();
     hinfo.dir.y = dir.Y();
     hinfo.dir.z = dir.Z();
 
     // And determine if the Hit is on a Calorimetry object
     for (const art::Ptr<anab::Calorimetry> &c: calo) {
       if (c->PlaneID().Plane != hinfo.h.plane) continue;
 
       // Found the plane! Now find the hit:
       for (unsigned i_calo = 0; i_calo < c->dQdx().size(); i_calo++) {
         if (c->TpIndices()[i_calo] == hkey) { // "TpIndices" match to the hit key
           // Fill the calo information associated with the hit
           hinfo.oncalo = true;
           hinfo.pitch = c->TrkPitchVec()[i_calo];
           hinfo.dqdx = c->dQdx()[i_calo];
           hinfo.rr = c->ResidualRange()[i_calo];
           break;
         }
       }
       break;
     }
   }
 
   // Save SpacePoint information
   if (sp) {
     hinfo.h.sp.x = sp->position().x();
     hinfo.h.sp.y = sp->position().y();
     hinfo.h.sp.z = sp->position().z();
 
     hinfo.h.hasSP = true;
   }
   else {
     hinfo.h.hasSP = false;
   }
 
   return hinfo;
 }

TrackCaloSkimmer& sbn::TrackCaloSkimmer::operator= ( TrackCaloSkimmer const & )

delete

TrackCaloSkimmer& sbn::TrackCaloSkimmer::operator= ( TrackCaloSkimmer && )

delete

void sbn::TrackCaloSkimmer::respondToOpenInputFile ( const art::FileBlock & fb )

inlineoverride

Definition at line 89 of file TrackCaloSkimmer.h.

                                                                {
     (void) fb;
     fMeta.ifile ++;
   }

Member Data Documentation

art::InputTag sbn::TrackCaloSkimmer::fCALOproducer

private

Definition at line 172 of file TrackCaloSkimmer.h.

bool sbn::TrackCaloSkimmer::fDoTailFit

private

Definition at line 180 of file TrackCaloSkimmer.h.

bool sbn::TrackCaloSkimmer::fFillTrackEndHits

private

Definition at line 186 of file TrackCaloSkimmer.h.

TFitter sbn::TrackCaloSkimmer::fFitConst

private

Definition at line 205 of file TrackCaloSkimmer.h.

TFitter sbn::TrackCaloSkimmer::fFitExp

private

Definition at line 204 of file TrackCaloSkimmer.h.

std::string sbn::TrackCaloSkimmer::fG4producer

private

Definition at line 177 of file TrackCaloSkimmer.h.

art::InputTag sbn::TrackCaloSkimmer::fHITproducer

private

Definition at line 175 of file TrackCaloSkimmer.h.

double sbn::TrackCaloSkimmer::fHitRawDigitsTickCollectWidth

private

Definition at line 183 of file TrackCaloSkimmer.h.

int sbn::TrackCaloSkimmer::fHitRawDigitsWireCollectWidth

private

Definition at line 185 of file TrackCaloSkimmer.h.

MetaInfo sbn::TrackCaloSkimmer::fMeta

private

Definition at line 194 of file TrackCaloSkimmer.h.

art::InputTag sbn::TrackCaloSkimmer::fPFPproducer

private

Definition at line 170 of file TrackCaloSkimmer.h.

std::vector<art::InputTag> sbn::TrackCaloSkimmer::fRawDigitproducers

private

Definition at line 176 of file TrackCaloSkimmer.h.

bool sbn::TrackCaloSkimmer::fRequireT0

private

Definition at line 179 of file TrackCaloSkimmer.h.

std::vector<std::unique_ptr<sbn::ITCSSelectionTool> > sbn::TrackCaloSkimmer::fSelectionTools

private

Definition at line 191 of file TrackCaloSkimmer.h.

bool sbn::TrackCaloSkimmer::fSilenceMissingDataProducts

private

Definition at line 182 of file TrackCaloSkimmer.h.

std::string sbn::TrackCaloSkimmer::fSimChannelproducer

private

Definition at line 178 of file TrackCaloSkimmer.h.

std::map<Snippet, int> sbn::TrackCaloSkimmer::fSnippetCount

private

Definition at line 196 of file TrackCaloSkimmer.h.

std::vector<art::InputTag> sbn::TrackCaloSkimmer::fT0producers

private

Definition at line 171 of file TrackCaloSkimmer.h.

double sbn::TrackCaloSkimmer::fTailFitResidualRange

private

Definition at line 184 of file TrackCaloSkimmer.h.

TrackInfo* sbn::TrackCaloSkimmer::fTrack

private

Definition at line 201 of file TrackCaloSkimmer.h.

float sbn::TrackCaloSkimmer::fTrackEndHitTimeBox

private

Definition at line 188 of file TrackCaloSkimmer.h.

float sbn::TrackCaloSkimmer::fTrackEndHitWireBox

private

Definition at line 187 of file TrackCaloSkimmer.h.

TTree* sbn::TrackCaloSkimmer::fTree

private

Definition at line 200 of file TrackCaloSkimmer.h.

art::InputTag sbn::TrackCaloSkimmer::fTRKHMproducer

private

Definition at line 174 of file TrackCaloSkimmer.h.

art::InputTag sbn::TrackCaloSkimmer::fTRKproducer

private

Definition at line 173 of file TrackCaloSkimmer.h.

bool sbn::TrackCaloSkimmer::fVerbose

private

Definition at line 181 of file TrackCaloSkimmer.h.

std::map<geo::WireID, std::pair<int, int> > sbn::TrackCaloSkimmer::fWiresToSave

private

Definition at line 197 of file TrackCaloSkimmer.h.

The documentation for this class was generated from the following files:

Classes

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation