Inheritance diagram for sbn::ShowerSelectionVars:

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

	ShowerSelectionVars (ShowerSelectionVars const &)=delete

	ShowerSelectionVars (ShowerSelectionVars &&)=delete

ShowerSelectionVars &	operator= (ShowerSelectionVars const &)=delete

ShowerSelectionVars &	operator= (ShowerSelectionVars &&)=delete

void	produce (art::Event &e) override

Private Member Functions
sbn::ShowerDensityFit	DensityFitter (const recob::Shower &shower, const std::vector< art::Ptr< recob::SpacePoint > > &sps) const

sbn::ShowerTrackFit	TrackFinder (const recob::Shower &shower, const recob::Track &track, const std::vector< art::Ptr< recob::SpacePoint > > &sps) const

std::vector< float >	ResidualFinder (const art::Ptr< recob::Shower > &primaryShower, const std::vector< art::Ptr< recob::Shower >> &sliceShowers) const

double	SpacePointPerpendicular (const recob::SpacePoint &sp, TVector3 const &vertex, TVector3 const &direction) const

double	SpacePointPerpendicular (const recob::SpacePoint &sp, TVector3 const &vertex, TVector3 const &direction, double const proj) const

double	SpacePointProjection (const recob::SpacePoint &sp, TVector3 const &vertex, TVector3 const &direction) const

TVector3	SpacePointPosition (recob::SpacePoint const &sp) const

Private Attributes
const art::InputTag	fPandoraLabel

const art::InputTag	fShowerLabel

const unsigned int	fNSegments

const bool	fRemoveStartFin

Detailed Description

Definition at line 40 of file ShowerSelectionVars_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 85 of file ShowerSelectionVars_module.cc.

   : EDProducer{p}
   , fPandoraLabel(p.get<art::InputTag>("PandoraLabel"))
   , fShowerLabel(p.get<art::InputTag>("ShowerLabel"))
   , fNSegments(p.get<unsigned int>("NSegments"))
   , fRemoveStartFin(p.get<bool>("RemoveStartFin"))
 {
   if (fNSegments==0)
     throw cet::exception("ShowerSelectionVars") << "Number of Segements must be non-zero" << std::endl;
 
   produces<std::vector<float> >();
   produces<std::vector<sbn::ShowerTrackFit> >();
   produces<std::vector<sbn::ShowerDensityFit> >();
 
   produces<art::Assns<recob::Shower, float> >();
   produces<art::Assns<recob::Shower, sbn::ShowerTrackFit> >();
   produces<art::Assns<recob::Shower, sbn::ShowerDensityFit> >();
 }

sbn::ShowerSelectionVars::ShowerSelectionVars ( ShowerSelectionVars const & )

delete

sbn::ShowerSelectionVars::ShowerSelectionVars ( ShowerSelectionVars && )

delete

Member Function Documentation

sbn::ShowerDensityFit sbn::ShowerSelectionVars::DensityFitter	(	const recob::Shower &	shower,
		const std::vector< art::Ptr< recob::SpacePoint > > &	sps
	)		const

private

Definition at line 234 of file ShowerSelectionVars_module.cc.

                                                           {
 
   const TVector3 showerVtx(shower.ShowerStart());
   const TVector3 showerDir(shower.Direction());
   const double OpenAngle(shower.OpenAngle());
   unsigned int totalHits(0);
 
   if (!shower.has_length() || !shower.has_open_angle() || sps.empty())
     return sbn::ShowerDensityFit();
 
   std::map<int, std::vector<art::Ptr<recob::SpacePoint> > > segmentSPMap;
   double segmentSize = shower.Length()/fNSegments;
 
   //Split the the spacepoints into segments.
   for(auto const& sp: sps){
 
     //Get the position of the spacepoint
     const double projLen(SpacePointProjection(*sp, showerVtx, showerDir));
     const double perpLen(SpacePointPerpendicular(*sp, showerVtx, showerDir, projLen));
 
     if(perpLen > TMath::Abs(TMath::Tan(OpenAngle)*projLen))
       continue;
 
     //Get where the sp should be place.
     const int sg_len(round(projLen/segmentSize));
     segmentSPMap[sg_len].push_back(sp);
     ++totalHits;
   }
 
   TGraph* graph = new TGraph();
 
   //Calculate the density gradent.
   for(auto& segment: segmentSPMap){
     double sg_len = segment.first;
 
     if(segment.second.size() < 10){continue;}
 
     //Calculate the charge in the segement
     double segmentHits = segment.second.size();
 
     //Calculate the voume
     double lower_dist = sg_len*segmentSize - segmentSize/2;
     double upper_dist = sg_len*segmentSize + segmentSize/2;
 
     if(fRemoveStartFin){if(sg_len==0 || sg_len==fNSegments){continue;}}
 
     if(sg_len==0)         {lower_dist = 0;}
     if(sg_len==fNSegments){upper_dist = sg_len*segmentSize;}
 
     double littlevolume = lower_dist*TMath::Power((TMath::Tan(0.5*OpenAngle)*lower_dist),2)*TMath::Pi()/3;
     double bigvolume    = upper_dist*TMath::Power((TMath::Tan(0.5*OpenAngle)*upper_dist),2)*TMath::Pi()/3;
     double volume       = bigvolume - littlevolume;
 
     double SegmentDensity = segmentHits/volume;
 
     double LengthToSegment = (lower_dist+upper_dist)/2;
 
     graph->SetPoint(graph->GetN(),LengthToSegment,SegmentDensity/totalHits);
   }
 
   if(graph->GetN() < 3 )
     return sbn::ShowerDensityFit();
 
   TF1 *fit = new TF1("fit", "[0]/x^[1]");
   fit->SetParLimits(1,1,2);
   fit->SetParLimits(0,0,1);
 
   graph->Fit(fit,"Q");
 
   const double grad(fit->GetParameter(0));
   const double pow(fit->GetParameter(1));
 
   delete fit;
   delete graph;
 
   return sbn::ShowerDensityFit(grad, pow);
 }

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

delete

ShowerSelectionVars& sbn::ShowerSelectionVars::operator= ( ShowerSelectionVars && )

delete

void sbn::ShowerSelectionVars::produce ( art::Event & e )

override

Definition at line 104 of file ShowerSelectionVars_module.cc.

 {
   //Get the showers
   auto const showerHandle = e.getValidHandle<std::vector<recob::Shower> >(fShowerLabel);
   auto const sliceHandle = e.getValidHandle<std::vector<recob::Slice> >(fPandoraLabel);
   auto const pfpHandle = e.getValidHandle<std::vector<recob::PFParticle> >(fPandoraLabel);
   // auto const spHandle = e.getValidHandle<std::vector<recob::SpacePoint> >(fPandoraLabel);
 
   std::vector<art::Ptr<recob::Shower> > showers;
   art::fill_ptr_vector(showers, showerHandle);
 
   std::vector<art::Ptr<recob::Slice> > slices;
   art::fill_ptr_vector(slices, sliceHandle);
 
   std::vector<art::Ptr<recob::PFParticle> > pfps;
   art::fill_ptr_vector(pfps, pfpHandle);
 
   art::FindManyP<recob::PFParticle> fmSlicePFP(sliceHandle, e, fPandoraLabel);
   if(!fmSlicePFP.isValid()){
     throw cet::exception("ShowerSelectionVars") << "Slice-PFP association is somehow not valid. Stopping";
     return;
   }
   art::FindManyP<recob::Shower> fmPFPShower(pfpHandle, e, fShowerLabel);
   if(!fmPFPShower.isValid()){
     throw cet::exception("ShowerSelectionVars") << "PFP-Shower association is somehow not valid. Stopping";
     return;
   }
   art::FindManyP<recob::SpacePoint> fmShowerSP(showerHandle, e, fShowerLabel);
   if(!fmShowerSP.isValid()){
     throw cet::exception("ShowerSelectionVars") << "Shower-SP association is somehow not valid. Stopping";
     return;
   }
   art::FindManyP<recob::SpacePoint> fmTrackSP(showerHandle, e, fShowerLabel);
   if(!fmTrackSP.isValid()){
     throw cet::exception("TrackSelectionVars") << "Track-SP association is somehow not valid. Stopping";
     return;
   }
   art::FindManyP<recob::Track> fmShowerTrack(showerHandle, e, fShowerLabel);
   if(!fmShowerTrack.isValid()){
     throw cet::exception("ShowerSelectionVars") << "Shower-Track association is somehow not valid. Stopping";
     return;
   }
   // art::FindManyP<recob::Track> fmSPHit(spHandle, e, fPandoraLabel);
   // if(!fmSPHit.isValid()){
   //   throw cet::exception("ShowerSelectionVars") << "SP-Hit association is somehow not valid. Stopping";
   //   return;
   // }
 
   // Create a map of PFP ID to PFP object for convinience
   std::map<size_t, art::Ptr<recob::PFParticle> > pfpMap;
   for (auto const& pfp: pfps) {
     pfpMap[pfp->Self()] = pfp;
   }
 
   std::map<art::Ptr<recob::Shower>, std::vector<float> > showerResidualMap;
 
   // Calculate the residuals: Want to calculate the residuals from each primary
   // shower to every other primary shower in the same slice
   for (auto const& slice: slices) {
     std::vector<art::Ptr<recob::PFParticle> >  slicePFPs(fmSlicePFP.at(slice.key()));
 
     if (slicePFPs.empty())
       continue;
 
     art::Ptr<recob::PFParticle> pfpNeutrino;
     for (auto const& pfp: slicePFPs) {
       if(std::abs(pfp->PdgCode()) == 12 || std::abs(pfp->PdgCode()) == 14){
         pfpNeutrino = pfp;
         break;
       } // if neutrino
     } // end pfp: slicePFPs
 
     if (pfpNeutrino.isNull())
       continue;
 
     std::vector<art::Ptr<recob::Shower>> primaryShowers;
     for (auto const& daughterID: pfpNeutrino->Daughters()) {
       const auto& pfp(pfpMap[daughterID]);
       if (pfp->PdgCode()==11) {
         primaryShowers.push_back(fmPFPShower.at(pfp.key()).front());
       } // if showers
     } // end daughterId: pfpNeutrino->Daughters
 
     if (primaryShowers.size()<2)
       continue;
 
     for (auto const& shower: primaryShowers){
       showerResidualMap[shower] = ResidualFinder(shower, primaryShowers);
     } // end shower: primaryShowers
   } // end slice: slices
 
   std::unique_ptr<std::vector<float> >residualCol(std::make_unique<std::vector<float> >());
   std::unique_ptr<std::vector<sbn::ShowerDensityFit> > densityFitCol(std::make_unique<std::vector<sbn::ShowerDensityFit> >());
   std::unique_ptr<std::vector<sbn::ShowerTrackFit> > trackFitCol(std::make_unique<std::vector<sbn::ShowerTrackFit> >());
 
   std::unique_ptr<art::Assns<recob::Shower, float> >
     residualAssns(std::make_unique<art::Assns<recob::Shower, float> >());
   std::unique_ptr<art::Assns<recob::Shower, sbn::ShowerDensityFit> >
     densityFitAssns(std::make_unique<art::Assns<recob::Shower, sbn::ShowerDensityFit> >());
   std::unique_ptr<art::Assns<recob::Shower, sbn::ShowerTrackFit> >
     trackFitAssns(std::make_unique<art::Assns<recob::Shower, sbn::ShowerTrackFit> >());
 
   for (auto const& shower: showers){
 
     for (auto const res: showerResidualMap[shower]) {
       residualCol->push_back(res);
       util::CreateAssn(*this, e, *residualCol, shower, *residualAssns);
     }
 
     sbn::ShowerDensityFit showerDensityFit(DensityFitter(*shower, fmShowerSP.at(shower.key())));
     densityFitCol->push_back(showerDensityFit);
     util::CreateAssn(*this, e, *densityFitCol, shower, *densityFitAssns);
 
     const std::vector<art::Ptr<recob::Track> > showerTrackVec(fmShowerTrack.at(shower.key()));
     sbn::ShowerTrackFit showerTrackFit(showerTrackVec.empty() ? sbn::ShowerTrackFit() :
         TrackFinder(*shower, *showerTrackVec.front(), fmTrackSP.at(showerTrackVec.front().key())));
     trackFitCol->push_back(showerTrackFit);
     util::CreateAssn(*this, e, *trackFitCol, shower, *trackFitAssns);
   }
 
   e.put(std::move(residualCol));
   e.put(std::move(densityFitCol));
   e.put(std::move(trackFitCol));
 
   e.put(std::move(residualAssns));
   e.put(std::move(densityFitAssns));
   e.put(std::move(trackFitAssns));
 }

std::vector< float > sbn::ShowerSelectionVars::ResidualFinder	(	const art::Ptr< recob::Shower > &	primaryShower,
		const std::vector< art::Ptr< recob::Shower >> &	sliceShowers
	)		const

private

Definition at line 335 of file ShowerSelectionVars_module.cc.

                                                                                                                                                           {
 
   std::vector<float> residuals;
 
   const TVector3 primaryVtx(primaryShower->ShowerStart());
   const TVector3 showerDir(primaryShower->Direction());
 
   // Check the vtx and direction are set
   if (primaryVtx.Z() < -500 || showerDir.Z() < -500)
     return residuals;
 
   for (auto const& shower: sliceShowers) {
 
     // Do not match shower to self
     if ((primaryVtx - shower->ShowerStart()).Mag() < std::numeric_limits<float>::epsilon())
       continue;
     // Check the vertex of the second shower is set
     if (shower->ShowerStart().Z() < -500) {
       residuals.push_back(-999);
       continue;
     }
 
     const TVector3 showerDisplacement(primaryVtx - shower->ShowerStart());
 
     const double projection(showerDir.Dot(showerDisplacement));
     const double perpendicular((showerDisplacement - projection*showerDir).Mag());
 
     residuals.push_back(perpendicular);
   }
 
   return residuals;
 }

double sbn::ShowerSelectionVars::SpacePointPerpendicular	(	const recob::SpacePoint &	sp,
		TVector3 const &	vertex,
		TVector3 const &	direction
	)		const

private

Definition at line 370 of file ShowerSelectionVars_module.cc.

                                                  {
 
   const double proj(SpacePointProjection(sp, vertex, direction));
   return SpacePointPerpendicular(sp, vertex, direction, proj);
  }

double sbn::ShowerSelectionVars::SpacePointPerpendicular	(	const recob::SpacePoint &	sp,
		TVector3 const &	vertex,
		TVector3 const &	direction,
		double const	proj
	)		const

private

Definition at line 378 of file ShowerSelectionVars_module.cc.

                                          {
 
   const TVector3 pos = SpacePointPosition(sp) - vertex;
   return (pos - proj*direction).Mag();
 }

TVector3 sbn::ShowerSelectionVars::SpacePointPosition ( recob::SpacePoint const & sp ) const

private

Definition at line 395 of file ShowerSelectionVars_module.cc.

                                                                                    {
 
   const Double32_t* sp_xyz = sp.XYZ();
   TVector3 sp_postiion = {sp_xyz[0], sp_xyz[1], sp_xyz[2]};
   return sp_postiion;
 }

double sbn::ShowerSelectionVars::SpacePointProjection	(	const recob::SpacePoint &	sp,
		TVector3 const &	vertex,
		TVector3 const &	direction
	)		const

private

Definition at line 387 of file ShowerSelectionVars_module.cc.

                                                  {
 
   TVector3 pos = SpacePointPosition(sp) - vertex;
   return pos.Dot(direction);
 }

sbn::ShowerTrackFit sbn::ShowerSelectionVars::TrackFinder	(	const recob::Shower &	shower,
		const recob::Track &	track,
		const std::vector< art::Ptr< recob::SpacePoint > > &	sps
	)		const

private

Definition at line 313 of file ShowerSelectionVars_module.cc.

                                                           {
 
   const unsigned int numSPs(sps.size());
 
   if (sps.empty())
     return sbn::ShowerTrackFit();
 
   double perp(0);
   const TVector3 showerVtx(shower.ShowerStart());
   const TVector3 showerDir(shower.Direction());
 
   for (auto const& sp: sps) {
     perp += SpacePointPerpendicular(*sp, showerVtx, showerDir);
   }
 
   const double trackLen(track.Length());
   const double trackWidth(perp/numSPs);
 
   return sbn::ShowerTrackFit(trackLen, trackWidth, numSPs);
 }