ana::StoppingParticleCosmicIdAlg Class Reference

#include <StoppingParticleCosmicIdAlg.h>

Classes
struct	Config

struct	Containment

Public Member Functions
	StoppingParticleCosmicIdAlg (const core::ProviderManager &manager, const Config &config)

	StoppingParticleCosmicIdAlg (const core::ProviderManager &manager, const fhicl::ParameterSet &pset)

	~StoppingParticleCosmicIdAlg ()

void	reconfigure (const core::ProviderManager &manager, const Config &config)

double	StoppingChiSq (geo::Point_t end, std::vector< art::Ptr< anab::Calorimetry >> calos)

bool	StoppingEnd (geo::Point_t end, std::vector< art::Ptr< anab::Calorimetry >> calos)

bool	StoppingParticleCosmicId (recob::Track track, std::vector< art::Ptr< anab::Calorimetry >> calos)

bool	StoppingParticleCosmicId (recob::Track track, recob::Track track2, std::vector< art::Ptr< anab::Calorimetry >> calos, std::vector< art::Ptr< anab::Calorimetry >> calos2)

bool	InFiducial (geo::Point_t point)

Private Attributes
double	fMinX

double	fMinY

double	fMinZ

double	fMaxX

double	fMaxY

double	fMaxZ

double	fResRangeMin

double	fResRangeMax

double	fDEdxMax

double	fStoppingChi2Limit

std::vector< geo::BoxBoundedGeo >	fFiducialVolumes

Detailed Description

Definition at line 36 of file sbnana/sbnanalysis/ana/SBNOscReco/CosmicIDAlgs/StoppingParticleCosmicIdAlg.h.

Constructor & Destructor Documentation

ana::StoppingParticleCosmicIdAlg::StoppingParticleCosmicIdAlg	(	const core::ProviderManager &	manager,
		const Config &	config
	)

Definition at line 5 of file sbnana/sbnanalysis/ana/SBNOscReco/CosmicIDAlgs/StoppingParticleCosmicIdAlg.cc.

                                                                                                                 {
 
   this->reconfigure(manager, config);
 
 }

ana::StoppingParticleCosmicIdAlg::StoppingParticleCosmicIdAlg	(	const core::ProviderManager &	manager,
		const fhicl::ParameterSet &	pset
	)

inline

Definition at line 85 of file sbnana/sbnanalysis/ana/SBNOscReco/CosmicIDAlgs/StoppingParticleCosmicIdAlg.h.

85 :

86 StoppingParticleCosmicIdAlg(manager, fhicl::Table<Config>(pset, {})()) {}

ana::StoppingParticleCosmicIdAlg::StoppingParticleCosmicIdAlg

StoppingParticleCosmicIdAlg(const core::ProviderManager &manager, const Config &config)

Definition: sbnana/sbnanalysis/ana/SBNOscReco/CosmicIDAlgs/StoppingParticleCosmicIdAlg.cc:5

ana::StoppingParticleCosmicIdAlg::~StoppingParticleCosmicIdAlg ( )

Definition at line 17 of file sbnana/sbnanalysis/ana/SBNOscReco/CosmicIDAlgs/StoppingParticleCosmicIdAlg.cc.

                                                          {
 
 }

Member Function Documentation

bool ana::StoppingParticleCosmicIdAlg::InFiducial ( geo::Point_t point )

Definition at line 151 of file sbnana/sbnanalysis/ana/SBNOscReco/CosmicIDAlgs/StoppingParticleCosmicIdAlg.cc.

                                                              {
   for (const geo::BoxBoundedGeo vol: fFiducialVolumes) {
     if (vol.ContainsPosition(point)) return true;
   }
   return false;
 }

void ana::StoppingParticleCosmicIdAlg::reconfigure	(	const core::ProviderManager &	manager,
		const Config &	config
	)

Definition at line 22 of file sbnana/sbnanalysis/ana/SBNOscReco/CosmicIDAlgs/StoppingParticleCosmicIdAlg.cc.

                                                                                                      {
   fFiducialVolumes.clear();
 
   fMinX = config.ContainmentCuts().MinX(); 
   fMinY = config.ContainmentCuts().MinY(); 
   fMinZ = config.ContainmentCuts().MinZ(); 
   fMaxX = config.ContainmentCuts().MaxX(); 
   fMaxY = config.ContainmentCuts().MaxY(); 
   fMaxZ = config.ContainmentCuts().MaxZ();
 
   // get the set of fiducial volumes
   for (auto const &cryo: manager.GetGeometryProvider()->IterateCryostats()) {
     geo::GeometryCore::TPC_iterator tstart = manager.GetGeometryProvider()->begin_TPC(cryo.ID()),
                                     tend = manager.GetGeometryProvider()->end_TPC(cryo.ID());
 
     double Min_X = std::min_element(tstart, tend, [](auto &lhs, auto &rhs) { return lhs.MinX() < rhs.MinX(); })->MinX();
     double Min_Y = std::min_element(tstart, tend, [](auto &lhs, auto &rhs) { return lhs.MinY() < rhs.MinY(); })->MinY();
     double Min_Z = std::min_element(tstart, tend, [](auto &lhs, auto &rhs) { return lhs.MinZ() < rhs.MinZ(); })->MinZ();
     
     double Max_X = std::max_element(tstart, tend, [](auto &lhs, auto &rhs) { return lhs.MaxX() < rhs.MaxX(); })->MaxX();
     double Max_Y = std::max_element(tstart, tend, [](auto &lhs, auto &rhs) { return lhs.MaxY() < rhs.MaxY(); })->MaxY();
     double Max_Z = std::max_element(tstart, tend, [](auto &lhs, auto &rhs) { return lhs.MaxZ() < rhs.MaxZ(); })->MaxZ();
 
     fFiducialVolumes.emplace_back(Min_X + fMinX, Max_X - fMaxX, Min_Y + fMinY, Max_Y - fMaxY, Min_Z + fMinZ, Max_Z - fMaxZ);
   }
 
   fResRangeMin = config.ResRangeMin();
   fResRangeMax = config.ResRangeMax();
   fDEdxMax = config.DEdxMax();
   fStoppingChi2Limit = config.StoppingChi2Limit();
 
   return;
 }

double ana::StoppingParticleCosmicIdAlg::StoppingChiSq	(	geo::Point_t	end,
		std::vector< art::Ptr< anab::Calorimetry >>	calos
	)

Definition at line 57 of file sbnana/sbnanalysis/ana/SBNOscReco/CosmicIDAlgs/StoppingParticleCosmicIdAlg.cc.

                                                                                                          {
 
   // If calorimetry object is null then return 0
   if(calos.size()==0) return -99999;
 
   // Loop over planes (Y->V->U) and choose the next plane's calorimetry if there are 1.5x more points (collection plane more reliable)
   size_t nhits = 0;
   art::Ptr<anab::Calorimetry> calo = calos[0];
   for( size_t i = calos.size(); i > 0; i--){
     if(calos[i-1]->dEdx().size() > nhits*1.5){
       nhits = calos[i-1]->dEdx().size();
       calo = calos[i-1];
     }
   }
 
   // If there's something wrong with the calorimetry object return null
   if(calo->XYZ().size() != nhits || nhits < 1) return -99999;
 
   // Get the distance from the track point and the start/end of calo data
   double distStart = (calo->XYZ()[0] - end).Mag2();
   double distEnd = (calo->XYZ()[nhits-1] - end).Mag2();
 
   double maxDedx = 0;
   double resrgStart = 0;
   std::vector<double> v_resrg;
   std::vector<double> v_dedx;
   // Loop over plane's calorimetry data
   for(size_t i = 0; i < nhits; i++){
     double dedx = calo->dEdx()[i];
     double resrg = calo->ResidualRange()[i];
 
     // Flip the residual range if the track and calo objects don't match up
     if(distStart < distEnd && calo->ResidualRange()[0] > calo->ResidualRange()[nhits-1]) resrg = calo->ResidualRange()[0] - calo->ResidualRange()[i];
     if(distStart > distEnd && calo->ResidualRange()[0] < calo->ResidualRange()[nhits-1]) resrg = calo->ResidualRange()[nhits-1] - calo->ResidualRange()[i];
 
     // Find the maximum dE/dx within a limit and the corresponding res range
     if(resrg < fResRangeMin && dedx > maxDedx && dedx < fDEdxMax){
       maxDedx = dedx;
       resrgStart = resrg;
     }
 
   }
 
   // Loop over it again
   for(size_t i = 0; i < nhits; i++){
     double dedx = calo->dEdx()[i];
     double resrg = calo->ResidualRange()[i];
 
     // Flip the residual range if the track and calo objects don't match up
     if(distStart < distEnd && calo->ResidualRange()[0] > calo->ResidualRange()[nhits-1]) resrg = calo->ResidualRange()[0] - calo->ResidualRange()[i];
     if(distStart > distEnd && calo->ResidualRange()[0] < calo->ResidualRange()[nhits-1]) resrg = calo->ResidualRange()[nhits-1] - calo->ResidualRange()[i];
 
     // Record all dE/dx and residual ranges below limits
     if(resrg > resrgStart && resrg < resrgStart + fResRangeMax && dedx < fDEdxMax){
       v_resrg.push_back(resrg);
       v_dedx.push_back(dedx);
     }
   }
 
   // Return null value if not enough points to do fits
   if(v_dedx.size() < 10) return -99999;
 
   // Try to do a pol0 fit
   TGraph *gdedx = new TGraph(v_dedx.size(), &v_resrg[0], &v_dedx[0]);
   try{ gdedx->Fit("pol0", "Q"); } catch(...){ delete gdedx; return -99999; }
   TF1* polfit = gdedx->GetFunction("pol0");
   double polchi2 = polfit->GetChisquare();
 
   // Try to do and exp fit
   try{ gdedx->Fit("expo", "Q"); } catch(...){ delete gdedx; return -99999; }
   TF1* expfit = gdedx->GetFunction("expo");
   double expchi2 = expfit->GetChisquare();
 
   delete gdedx;
 
   // Return the chi2 ratio
   return polchi2/expchi2;
 
 }

bool ana::StoppingParticleCosmicIdAlg::StoppingEnd	(	geo::Point_t	end,
		std::vector< art::Ptr< anab::Calorimetry >>	calos
	)

Definition at line 139 of file sbnana/sbnanalysis/ana/SBNOscReco/CosmicIDAlgs/StoppingParticleCosmicIdAlg.cc.

                                                                                                      {
   
   // Get the chi2 ratio
   double chiSqRatio = StoppingChiSq(end, calos);
 
   // If it is below a limit then tag it as stopping
   if(chiSqRatio > fStoppingChi2Limit) return true;
 
   return false;
 }

bool ana::StoppingParticleCosmicIdAlg::StoppingParticleCosmicId	(	recob::Track	track,
		std::vector< art::Ptr< anab::Calorimetry >>	calos
	)

Definition at line 159 of file sbnana/sbnanalysis/ana/SBNOscReco/CosmicIDAlgs/StoppingParticleCosmicIdAlg.cc.

                                                                                                                     {
 
   // Check if start and end of track is inside the fiducial volume
   bool startInFiducial = InFiducial(track.Vertex());
   bool endInFiducial = InFiducial(track.End());
 
   // Check if the start and end of track stops
   bool startStops = StoppingEnd(track.Vertex(), calos);
   bool endStops = StoppingEnd(track.End(), calos);
 
   // If one end stops and the other end is outside of the FV then tag as stopping cosmic if the track is going downwards
   if((startStops && !endInFiducial && track.End().Y() > track.Vertex().Y()) 
       || (endStops && !startInFiducial && track.Vertex().Y() > track.End().Y())){
     return true;
   }
 
   return false;
 
 }

bool ana::StoppingParticleCosmicIdAlg::StoppingParticleCosmicId	(	recob::Track	track,
		recob::Track	track2,
		std::vector< art::Ptr< anab::Calorimetry >>	calos,
		std::vector< art::Ptr< anab::Calorimetry >>	calos2
	)