Inheritance diagram for ShowerRecoTools::ShowerBayesianTrucatingdEdx:

Public Member Functions
	ShowerBayesianTrucatingdEdx (const fhicl::ParameterSet &pset)

int	CalculateElement (const art::Ptr< recob::PFParticle > &pfparticle, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder) override

Public Member Functions inherited from ShowerRecoTools::IShowerTool
	IShowerTool (const fhicl::ParameterSet &pset)

virtual	~IShowerTool () noexcept=default

int	RunShowerTool (const art::Ptr< recob::PFParticle > &pfparticle, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder, std::string evd_display_name_append="")

virtual void	InitialiseProducers ()

void	SetPtr (art::ProducesCollector *collector)

void	InitaliseProducerPtr (reco::shower::ShowerProducedPtrsHolder &uniqueproducerPtrs)

virtual int	AddAssociations (const art::Ptr< recob::PFParticle > &pfpPtr, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder)

Private Member Functions
double	CalculatePosterior (std::string priorname, std::vector< double > &values, int &minprob, float &mean, float &likelihood)

double	CalculatePosterior (std::string priorname, std::vector< double > &values)

bool	isProbabilityGood (float &old_prob, float &new_prob)

bool	isPosteriorProbabilityGood (double &prob, double &old_posteior)

bool	CheckPoint (std::string priorname, double &value)

std::vector< double >	GetLikelihooddEdxVec (double &electronprob, double &photonprob, std::string prior, std::vector< double > &dEdxVec)

std::vector< double >	MakeSeed (std::vector< double > &dEdxVec)

void	ForceSeedToFit (std::vector< double > &SeedTrack, std::string &prior, float &mean, double &posterior)

void	RecurivelyAddHit (std::vector< double > &SeedTrack, std::vector< double > &dEdxVec, std::string &prior, int &SkippedHitsNum, float &old_mean, double &old_posteior)

Private Attributes
TH1F *	electronpriorHist

TH1F *	photonpriorHist

int	fVerbose

std::string	fdEdxInputLabel

int	fNumSeedHits

float	fProbSeedCut

float	fProbPointCut

float	fPostiorCut

int	fnSkipHits

std::string	fShowerdEdxOutputLabel

bool	fDefineBestPlane

std::string	fShowerBestPlaneOutputLabel

Additional Inherited Members
Protected Member Functions inherited from ShowerRecoTools::IShowerTool
const shower::LArPandoraShowerAlg &	GetLArPandoraShowerAlg () const

template<class T >
art::Ptr< T >	GetProducedElementPtr (std::string Name, reco::shower::ShowerElementHolder &ShowerEleHolder, int iter=-1)

template<class T >
void	InitialiseProduct (std::string Name, std::string InstanceName="")

template<class T , class A , class B >
void	AddSingle (A &a, B &b, std::string Name)

int	GetVectorPtrSize (std::string Name)

void	PrintPtrs ()

void	PrintPtr (std::string Name)

Detailed Description

Definition at line 26 of file ShowerBayesianTrucatingdEdx_tool.cc.

Constructor & Destructor Documentation

ShowerRecoTools::ShowerBayesianTrucatingdEdx::ShowerBayesianTrucatingdEdx ( const fhicl::ParameterSet & pset )

Definition at line 93 of file ShowerBayesianTrucatingdEdx_tool.cc.

     : IShowerTool(pset.get<fhicl::ParameterSet>("BaseTools"))
     , fVerbose(pset.get<int>("Verbose"))
     , fdEdxInputLabel(pset.get<std::string>("dEdxInputLabel"))
     , fNumSeedHits(pset.get<int>("NumSeedHits"))
     , fProbSeedCut(pset.get<float>("ProbSeedCut"))
     , fProbPointCut(pset.get<float>("ProbPointCut"))
     , fPostiorCut(pset.get<float>("PostiorCut"))
     , fnSkipHits(pset.get<int>("nSkipHits"))
     , fShowerdEdxOutputLabel(pset.get<std::string>("ShowerdEdxOutputLabel"))
     , fDefineBestPlane(pset.get<bool>("DefineBestPlane"))
     , fShowerBestPlaneOutputLabel(pset.get<std::string>("ShowerBestPlaneOutputLabel"))
   {
 
     //Get the prior file name
     std::string fname;
     cet::search_path sp("FW_SEARCH_PATH");
     auto PriorPath = pset.get<std::string>("PriorFname");
     if (!sp.find_file(PriorPath, fname)) {
       throw cet::exception("ShowerBayesianTrucatingdEdx") << "Could not find the prior file";
     }
     std::string electron_histoname = pset.get<std::string>("PriorElectronHistoName");
     std::string photon_histoname = pset.get<std::string>("PriorPhotonHistoName");
 
     TFile fin(fname.c_str(), "READ");
     if (!fin.IsOpen()) {
       throw cet::exception("ShowerBayesianTrucatingdEdx") << "Could read the prior file. Stopping";
     }
 
     //Get the histograms.
     electronpriorHist = dynamic_cast<TH1F*>(fin.Get(electron_histoname.c_str()));
     if (!electronpriorHist) {
       throw cet::exception("ShowerBayesianTrucatingdEdx") << "Could not read the electron hist";
     }
     photonpriorHist = dynamic_cast<TH1F*>(fin.Get(photon_histoname.c_str()));
     if (!photonpriorHist) {
       throw cet::exception("ShowerBayesianTrucatingdEdx") << "Could not read the photon hist ";
     }
 
     if (electronpriorHist->GetNbinsX() != photonpriorHist->GetNbinsX()) {
       throw cet::exception("ShowerBayesianTrucatingdEdx") << "Histrogram bins do not match";
     }
 
     //Normalise the histograms.
     electronpriorHist->Scale(1 / electronpriorHist->Integral());
     photonpriorHist->Scale(1 / photonpriorHist->Integral());
   }

Member Function Documentation

int ShowerRecoTools::ShowerBayesianTrucatingdEdx::CalculateElement	(	const art::Ptr< recob::PFParticle > &	pfparticle,
		art::Event &	Event,
		reco::shower::ShowerElementHolder &	ShowerEleHolder
	)

overridevirtual

Implements ShowerRecoTools::IShowerTool.

Definition at line 142 of file ShowerBayesianTrucatingdEdx_tool.cc.

   {
 
     //The idea , to some peoples distaste, is to attempt to improve the dEdx value by assuming
     //The particle is either a) electron b) a e-e+ pair.
     //We will take the start of track and work down until a few hits destory our postier probability.
 
     //Note: tried takeing the postior with the highest sum of the probabilitys on all three
     //      planes and on the 2 planes with the most hits. Made things worse.
 
     //Get the vectors of the dEdx Elements
     if (!ShowerEleHolder.CheckElement(fdEdxInputLabel)) {
       if (fVerbose)
         mf::LogError("ShowerBayesianTrucatingdEdx")
           << "Start position not set, returning " << std::endl;
       return 1;
     }
 
     std::map<int, std::vector<double>> dEdx_plane_final;
     std::map<int, std::vector<double>> dEdx_vec_planes;
     ShowerEleHolder.GetElement(fdEdxInputLabel, dEdx_vec_planes);
 
     //Do this for each plane;
     for (auto const& dEdx_vec_plane : dEdx_vec_planes) {
 
       //Set up out final value if we don't have any points.
       if (dEdx_vec_plane.second.size() < 1) {
         dEdx_plane_final[dEdx_vec_plane.first] = {};
         continue;
       }
 
       std::vector<double> dEdx_vec = dEdx_vec_plane.second;
 
       double electronprob_eprior = 0;
       double photonprob_eprior = 0;
 
       double electronprob_pprior = 0;
       double photonprob_pprior = 0;
 
       std::vector<double> dEdx_electronprior =
         GetLikelihooddEdxVec(electronprob_eprior, photonprob_eprior, "electron", dEdx_vec);
       std::vector<double> dEdx_photonprior =
         GetLikelihooddEdxVec(electronprob_pprior, photonprob_pprior, "photon", dEdx_vec);
 
       //Use the vector which maximises both priors.
       if (electronprob_eprior < photonprob_pprior) {
         dEdx_plane_final[dEdx_vec_plane.first] = dEdx_photonprior;
       }
       else {
         dEdx_plane_final[dEdx_vec_plane.first] = dEdx_electronprior;
       }
 
     } //Plane Loop
 
     //Calculate the median of the of dEdx.
     std::vector<double> dEdx_final;
     std::vector<double> dEdx_finalErr;
 
     int max_hits = -999;
     int best_plane = -999;
 
     bool check = false;
     for (auto const& dEdx_plane : dEdx_plane_final) {
 
       //Redefine the best plane
       if ((int)(dEdx_plane.second).size() > max_hits) {
         best_plane = dEdx_plane.first;
         max_hits = (dEdx_plane.second).size();
       }
 
       if ((dEdx_plane.second).empty()) {
         dEdx_final.push_back(-999);
         dEdx_finalErr.push_back(-999);
         continue;
       }
 
       dEdx_final.push_back(TMath::Median((dEdx_plane.second).size(), &(dEdx_plane.second)[0]));
       dEdx_finalErr.push_back(-999);
       check = true;
     }
 
     //Check at least one plane has the information
     if (!check) return 1;
 
     if (fDefineBestPlane) { ShowerEleHolder.SetElement(best_plane, fShowerBestPlaneOutputLabel); }
 
     ShowerEleHolder.SetElement(dEdx_final, dEdx_finalErr, fShowerdEdxOutputLabel);
 
     return 0;
   }

double ShowerRecoTools::ShowerBayesianTrucatingdEdx::CalculatePosterior	(	std::string	priorname,
		std::vector< double > &	values,
		int &	minprob,
		float &	mean,
		float &	likelihood
	)

private

Definition at line 246 of file ShowerBayesianTrucatingdEdx_tool.cc.

   {
 
     //Posterior prob;
     float posterior = 1;
     float meanprob = 0;
     float likelihood_other = 1;
     likelihood = 1;
 
     //Minimum probability temp
     float minprob_temp = 9999;
     minprob_iter = 0;
 
     TH1F* prior_hist = nullptr;
     TH1F* other_hist = nullptr;
 
     if (priorname == "electron") {
       prior_hist = electronpriorHist;
       other_hist = photonpriorHist;
     }
     if (priorname == "photon") {
       prior_hist = photonpriorHist;
       other_hist = electronpriorHist;
     }
 
     TAxis* xaxis = prior_hist->GetXaxis();
 
     //Loop over the hits and calculate the probability
     for (int i = 0; i < (int)values.size(); ++i) {
 
       float value = values[i];
 
       Int_t bin = xaxis->FindBin(value);
 
       float prob = -9999;
       float other_prob = -9999;
 
       if (bin != xaxis->GetNbins() || bin == 0) {
         //Calculate the likelihood
         prob = prior_hist->GetBinContent(bin);
         other_prob = other_hist->GetBinContent(bin);
       }
       else {
         prob = 0;
         other_prob = 0;
       }
 
       if (prob < minprob_temp) {
         minprob_temp = prob;
         minprob_iter = i;
       }
 
       if (prob == 0 && other_prob == 0) { continue; }
 
       //Calculate the posterior the mean probability and liklihood
       meanprob += prior_hist->GetBinContent(bin);
       likelihood *= prob;
       likelihood_other *= other_prob;
     }
 
     posterior = likelihood / (likelihood + likelihood_other);
 
     meanprob /= values.size();
     mean = meanprob;
     return posterior;
   }

double ShowerRecoTools::ShowerBayesianTrucatingdEdx::CalculatePosterior	(	std::string	priorname,
		std::vector< double > &	values
	)

private

Definition at line 236 of file ShowerBayesianTrucatingdEdx_tool.cc.

   {
     int minprob_iter = -999;
     float mean = -999;
     float likelihood = -999;
     return CalculatePosterior(priorname, values, minprob_iter, mean, likelihood);
   }

bool ShowerRecoTools::ShowerBayesianTrucatingdEdx::CheckPoint	(	std::string	priorname,
		double &	value
	)

private

Definition at line 318 of file ShowerBayesianTrucatingdEdx_tool.cc.

   {
 
     TH1F* prior_hist = nullptr;
 
     if (priorname == "electron") { prior_hist = electronpriorHist; }
     if (priorname == "photon") { prior_hist = photonpriorHist; }
 
     TAxis* xaxis = prior_hist->GetXaxis();
 
     Int_t bin = xaxis->FindBin(value);
 
     float prob = -9999;
 
     if (bin != xaxis->GetNbins() + 1 || bin == 0) {
       //Calculate the likelihood
       prob = prior_hist->GetBinContent(bin);
     }
     else {
       prob = 0;
     }
 
     //Return the probability of getting that point.
     return prob > fProbPointCut;
   }

void ShowerRecoTools::ShowerBayesianTrucatingdEdx::ForceSeedToFit	(	std::vector< double > &	SeedTrack,
		std::string &	prior,
		float &	mean,
		double &	posterior
	)

private

Definition at line 399 of file ShowerBayesianTrucatingdEdx_tool.cc.

   {
 
     int minprob_iter = 999;
     float likelihood = -999;
     float prob = CalculatePosterior(prior, SeedTrack, minprob_iter, mean, likelihood);
     while ((mean < fProbSeedCut || prob <= 0) && SeedTrack.size() > 1) {
 
       //Remove the the worse point.
       // std::cout << "removing hit with dEdx: " << SeedTrack.at(minprob_iter) << std::endl;
       SeedTrack.erase(SeedTrack.begin() + minprob_iter);
       minprob_iter = 999;
 
       //Recalculate
       prob = CalculatePosterior(prior, SeedTrack, minprob_iter, mean, likelihood);
     }
     posterior = prob;
     // std::cout << "seed has been fit with size: " << SeedTrack.size() << std::endl;
     return;
   }

std::vector< double > ShowerRecoTools::ShowerBayesianTrucatingdEdx::GetLikelihooddEdxVec	(	double &	electronprob,
		double &	photonprob,
		std::string	prior,
		std::vector< double > &	dEdxVec
	)

private

Definition at line 345 of file ShowerBayesianTrucatingdEdx_tool.cc.

   {
 
     //have a pool
     std::vector<double> dEdxVec_temp = dEdxVec;
 
     //Get The seed track.
     std::vector<double> SeedTrack = MakeSeed(dEdxVec_temp);
     //    if(SeedTrack.size() < 1){
     //  return SeedTrack;
     // }
 
     //Force the seed the be a good likelihood.
     float mean = 999;
     double posterior = 999;
     ForceSeedToFit(SeedTrack, prior, mean, posterior);
 
     //Recursively add dEdx
     int SkippedHitsNum = 0;
     RecurivelyAddHit(SeedTrack, dEdxVec_temp, prior, SkippedHitsNum, mean, posterior);
 
     //Calculate the likelihood of the vector  with the photon and electron priors.
     electronprob = CalculatePosterior("electron", SeedTrack);
     photonprob = CalculatePosterior("photon", SeedTrack);
 
     return SeedTrack;
   }

bool ShowerRecoTools::ShowerBayesianTrucatingdEdx::isPosteriorProbabilityGood	(	double &	prob,
		double &	old_posteior
	)

inlineprivate

Definition at line 51 of file ShowerBayesianTrucatingdEdx_tool.cc.

     {
       return (old_posteior - prob) < fPostiorCut;
     }

bool ShowerRecoTools::ShowerBayesianTrucatingdEdx::isProbabilityGood	(	float &	old_prob,
		float &	new_prob
	)

inlineprivate

Definition at line 45 of file ShowerBayesianTrucatingdEdx_tool.cc.

     {
       return (old_prob - new_prob) < fProbSeedCut;
     }

std::vector< double > ShowerRecoTools::ShowerBayesianTrucatingdEdx::MakeSeed ( std::vector< double > & dEdxVec )

private

Definition at line 377 of file ShowerBayesianTrucatingdEdx_tool.cc.

   {
 
     std::vector<double> seed_vector;
 
     //Add the first hits to the seed
     int MaxHit = fNumSeedHits;
     if (fNumSeedHits > (int)dEdxVec.size()) { MaxHit = (int)dEdxVec.size(); }
 
     //    if(MaxHit == 0){
     //  mf::LogError("ShowerBayesianTrucatingdEdx") << "Size of the vector is 0 cannot perform the dEdx cutting "<< std::endl;
     //}
 
     for (int hit_iter = 0; hit_iter < MaxHit; ++hit_iter) {
       seed_vector.push_back(dEdxVec[0]);
       dEdxVec.erase(dEdxVec.begin());
     }
 
     return seed_vector;
   }

void ShowerRecoTools::ShowerBayesianTrucatingdEdx::RecurivelyAddHit	(	std::vector< double > &	SeedTrack,
		std::vector< double > &	dEdxVec,
		std::string &	prior,
		int &	SkippedHitsNum,
		float &	old_mean,
		double &	old_posteior
	)

private

Definition at line 424 of file ShowerBayesianTrucatingdEdx_tool.cc.

   {
 
     //If we have no more hits to add then lets finish.
     if (dEdxVec.size() < 1) { return; }
 
     bool ok = CheckPoint(prior, dEdxVec[0]);
     // int minprob_iter=999;
     // float mean = -999;
     // float likelihood =999;
     // if(ok){std::cout << "passed first cut" << std::endl;}
     // else{std::cout << "failed first cut" << std::endl;}
     // double prob = CalculatePosterior(prior,SeedTrack,minprob_iter,mean,likelihood);
     // ok *= isProbabilityGood(mean,old_mean);
     // if(ok){std::cout << "passed second cut" << std::endl;}
     // else{std::cout << "failed second cut" << std::endl;}
     // ok *= isPosteriorProbabilityGood(prob,old_posteior);
     // if(ok){std::cout << "passed this cut" << std::endl;}
     // else{std::cout << "failed third cut" << std::endl;}
 
     //If we failed lets try the next hits
     if (!ok) {
       // std::cout << "failed the pass point: " << dEdxVec[0] << " trying another hit" << SkippedHitsNum << std::endl;
       //if(SeedTrack.size() > 1){SeedTrack.pop_back();}
       ++SkippedHitsNum;
       if (SkippedHitsNum > fnSkipHits) { return; }
     }
     else {
       //Add the next point in question.
       // std::cout << "adding value: " << dEdxVec[0] << std::endl;
       //Reset the skip number
       SkippedHitsNum = 0;
       SeedTrack.push_back(dEdxVec[0]);
     }
 
     //We have analysed this hit now lets remove it.
     dEdxVec.erase(dEdxVec.begin());
 
     RecurivelyAddHit(SeedTrack, dEdxVec, prior, SkippedHitsNum, old_mean, old_posteior);
 
     return;
   }