Inheritance diagram for sbn::PCAngleKinkFinder:

Classes
struct	FitResult

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

	PCAngleKinkFinder (PCAngleKinkFinder const &)=delete

	PCAngleKinkFinder (PCAngleKinkFinder &&)=delete

PCAngleKinkFinder &	operator= (PCAngleKinkFinder const &)=delete

PCAngleKinkFinder &	operator= (PCAngleKinkFinder &&)=delete

void	produce (art::Event &e) override

Private Member Functions
FitResult	Fit (const std::vector< sbn::PCAngle > &angles, unsigned lo, unsigned hi, unsigned center, float est_angle, float est_pitch)

Private Attributes
std::string	fPFParticleLabel

std::string	fAngleLabel

float	fHitThreshold

float	fMinDist

bool	fAllowIncomplete

float	fPCARadius

float	fFitResolution

Detailed Description

Definition at line 46 of file PCAngleKinkFinder_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 275 of file PCAngleKinkFinder_module.cc.

   : EDProducer{p},
     fPFParticleLabel(p.get<std::string>("PFParticleLabel")),
     fAngleLabel(p.get<std::string>("AngleLabel")),
     fHitThreshold(p.get<float>("HitThreshold")),
     fMinDist(p.get<float>("MinDist")),
     fAllowIncomplete(p.get<bool>("AllowIncomplete")),
     fPCARadius(p.get<float>("PCARadius")),
     fFitResolution(p.get<float>("FitResolution", 0.05))
 {
   produces<std::vector<sbn::PCAngleKink>>();
   produces<art::Assns<recob::PFParticle, sbn::PCAngleKink>>();
 
 }

sbn::PCAngleKinkFinder::PCAngleKinkFinder ( PCAngleKinkFinder const & )

delete

sbn::PCAngleKinkFinder::PCAngleKinkFinder ( PCAngleKinkFinder && )

delete

Member Function Documentation

sbn::PCAngleKinkFinder::FitResult sbn::PCAngleKinkFinder::Fit	(	const std::vector< sbn::PCAngle > &	angles,
		unsigned	lo,
		unsigned	hi,
		unsigned	center,
		float	est_angle,
		float	est_pitch
	)

private

Definition at line 331 of file PCAngleKinkFinder_module.cc.

                                                                                                                                                                         {
   // setup the fitdata
   fFitDataX.clear();
   fFitDataY.clear();
   for (unsigned i = lo; i <= hi; i++) {
     fFitDataX.push_back( sqrt((angles[i].hit_time_cm - angles[center].hit_time_cm) * (angles[i].hit_time_cm - angles[center].hit_time_cm) +
       (angles[i].hit_wire_cm - angles[center].hit_wire_cm) * (angles[i].hit_wire_cm - angles[center].hit_wire_cm)) );
     fFitDataY.push_back(angles[i].angle);
   }
   fFitCenter = center;
 
   // setup the fitter
   TFitter fitter(5); // 5 param
   fitter.SetFCN(&FitLikelihood);
   fitter.SetParameter(0, "angle", est_angle, 0.1, 0, M_PI);
   fitter.SetParameter(1, "resolution", fFitResolution, 0., fFitResolution, fFitResolution);
   fitter.SetParameter(2, "radius", fPCARadius, 0., fPCARadius, fPCARadius);
   fitter.SetParameter(3, "pitch", est_pitch, 0., est_pitch, est_pitch);
   fitter.GetMinuit()->FixParameter(1);
   fitter.GetMinuit()->FixParameter(2);
   fitter.GetMinuit()->FixParameter(3);
 
   fitter.ExecuteCommand("MIGRAD", 0, 0);
 
   sbn::PCAngleKinkFinder::FitResult ret;
 
   ret.angle = fitter.GetParameter(0);
   ret.pitch = fitter.GetParameter(3);
 
   double chi2;
   double edm;
   double errdef;
   int nvpar;
   int nparx;
 
   fitter.GetStats(chi2, edm, errdef, nvpar, nparx);
   std::cout << "Fit chi2: " << chi2 << std::endl;
   std::cout << "Fit angle: " << ret.angle << std::endl;
 
   ret.chi2 = chi2;
 
   return ret;
 }

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

delete

PCAngleKinkFinder& sbn::PCAngleKinkFinder::operator= ( PCAngleKinkFinder && )

delete

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

override

Definition at line 375 of file PCAngleKinkFinder_module.cc.

 {
   // output stuff
   std::unique_ptr<std::vector<sbn::PCAngleKink>> outKinks(new std::vector<sbn::PCAngleKink>);
   std::unique_ptr<art::Assns<recob::PFParticle, sbn::PCAngleKink>> assn(new art::Assns<recob::PFParticle, sbn::PCAngleKink>);
 
   art::PtrMaker<sbn::PCAngleKink> kinkPtrMaker {evt};
 
   // get the PFParticle's and the associated data
   art::Handle<std::vector<recob::PFParticle>> pfparticle_handle;
   evt.getByLabel(fPFParticleLabel, pfparticle_handle);
 
   std::vector<art::Ptr<recob::PFParticle>> pfparticles;
   art::fill_ptr_vector(pfparticles, pfparticle_handle);
 
   // organize the PFPlist into a map
   std::map<unsigned, art::Ptr<recob::PFParticle>> id_to_pfp;
   for (unsigned i = 0; i < pfparticles.size(); i++) {
     id_to_pfp[pfparticles[i]->Self()] = pfparticles[i];
   }
 
   art::FindManyP<sbn::PCAnglePlane> pfparticleAngles(pfparticles, evt, fAngleLabel);
 
   // iterate over each particle
   for (unsigned i_part = 0; i_part < pfparticles.size(); i_part++) {
     art::Ptr<recob::PFParticle> thisPFP = pfparticles[i_part];
     std::cout << "NEW Particle: " << thisPFP->Self() << std::endl;
     const std::vector<art::Ptr<sbn::PCAnglePlane>> &thisAngles = pfparticleAngles.at(i_part);
 
     if (!thisAngles.size()) continue;
 
     std::vector<sbn::PCAngleKink> thisKinks;
     for (const art::Ptr<sbn::PCAnglePlane> &angle: thisAngles) {
       for (unsigned i_branch = 0; i_branch < angle->nBranches; i_branch++) {
         // flatten the structure of this branch
         auto [branchAngles, branch_start] = FlattenBranch(*angle, i_branch, fAllowIncomplete, fMinDist);
         bool in_hit = false;
         int hit_start_ind = -1;
         for (unsigned i_hit = branch_start; i_hit < branchAngles.size(); i_hit++) {
           bool this_in_hit = branchAngles[i_hit].angle > fHitThreshold;
           // New hit! set it up
           if (!in_hit && this_in_hit) {
             hit_start_ind = i_hit;
           }
           // end of hit! finish it up
           else if (in_hit && (!this_in_hit || i_hit+1 == branchAngles.size())) {
             float height = EstimateHeight(branchAngles, hit_start_ind, i_hit);
 
             int hi_hit_ind = FindFWHMHiIndex(branchAngles, hit_start_ind, i_hit, height, fHitThreshold);
             int lo_hit_ind = FindFWHMLoIndex(branchAngles, hit_start_ind, i_hit, height, fHitThreshold);
             // invalid hit
             // check if hit is valid
             if (hi_hit_ind >= 0 && lo_hit_ind >= 0) {
               std::cout << "NEW KINK: Particle: " << thisPFP->Self() << " Branch Ind: " << i_branch << " Plane: " << angle->plane.Plane << std::endl;
               std::cout << "EST HEIGHT: " << height << " IND LO: " << lo_hit_ind << " IND HI: " << hi_hit_ind << std::endl;
 
               // build out the kink information
               int max_ind = FindMaxIndex(branchAngles, lo_hit_ind, hi_hit_ind);
 
               // estimate the track pitch -- distance between successive hits
               float dist_lo_sum = 0.;
               for (int i = lo_hit_ind; i < max_ind; i++) {
                 dist_lo_sum += sqrt((branchAngles[i].hit_time_cm - branchAngles[i+1].hit_time_cm) * (branchAngles[i].hit_time_cm - branchAngles[i+1].hit_time_cm) +
                     (branchAngles[i].hit_wire_cm - branchAngles[i+1].hit_wire_cm) * (branchAngles[i].hit_wire_cm - branchAngles[i+1].hit_wire_cm)); 
               }
               // float pitch_lo = dist_lo_sum / (max_ind - lo_hit_ind);
 
               float dist_hi_sum = 0.;
               for (int i = max_ind; i < hi_hit_ind; i++) {
                 dist_hi_sum += sqrt((branchAngles[i].hit_time_cm - branchAngles[i+1].hit_time_cm) * (branchAngles[i].hit_time_cm - branchAngles[i+1].hit_time_cm) +
                     (branchAngles[i].hit_wire_cm - branchAngles[i+1].hit_wire_cm) * (branchAngles[i].hit_wire_cm - branchAngles[i+1].hit_wire_cm)); 
               }
               // float pitch_hi = dist_hi_sum / (hi_hit_ind - max_ind);
               float pitch = (dist_lo_sum + dist_hi_sum) / (hi_hit_ind - lo_hit_ind);
 
 
               sbn::PCAngleKinkFinder::FitResult fit = Fit(branchAngles, lo_hit_ind, hi_hit_ind, max_ind, height, pitch);
 
               const sbn::PCAngle &max = branchAngles[max_ind];
               const sbn::PCAngle &lo = branchAngles[lo_hit_ind];
               const sbn::PCAngle &hi = branchAngles[hi_hit_ind];
 
               sbn::PCAngleKink kink = BuildKink(max, lo, hi, fit, height);
               thisKinks.push_back(kink);
               // update the hit index to the end of this hit
               i_hit = hi_hit_ind;
             }
           }
 
           in_hit = this_in_hit;
         } 
       }
     }
 
     // save the kinks
     for (const sbn::PCAngleKink &k: thisKinks) {
       outKinks->push_back(k);
       art::Ptr<sbn::PCAngleKink> thisKinkPtr = kinkPtrMaker(outKinks->size()-1);
       assn->addSingle(thisPFP, thisKinkPtr);
     }
   }
 
   evt.put(std::move(outKinks));
   evt.put(std::move(assn));
 
 }