pid::PIDAAlg Class Reference

#include <PIDAAlg.h>

Classes
struct	PIDAProperties

Public Member Functions
	PIDAAlg (fhicl::ParameterSet const &p)
void	RunPIDAAlg (std::vector< float > const &, std::vector< float > const &)
void	RunPIDAAlg (anab::Calorimetry const &)
void	RunPIDAAlg (anab::Calorimetry const &, float &, float &)
float	getPIDAMean ()
float	getPIDASigma ()
const std::vector< float > &	getPIDAValues ()
const std::vector< float > &	getPIDAErrors ()
size_t	getNKDEBandwidths ()
float	getKDEBandwidth (const size_t i_b)
float	getPIDAKDEMostProbable (const size_t)
float	getPIDAKDEFullWidthHalfMax (const size_t)
void	PrintPIDAValues ()
void	setExponentConstant (float const &ex)
void	SetPIDATree (TTree *, TH1F *, std::vector< TH1F * >)
void	FillPIDATree (unsigned int, unsigned int, unsigned int, anab::Calorimetry const &)

Private Types
typedef struct pid::PIDAAlg::PIDAProperties	PIDAProperties_t

Private Member Functions
void	calculatePIDAMean ()
void	calculatePIDASigma ()
void	calculatePIDAIntegral (std::map< double, double > const &)
void	ClearInternalData ()
void	createKDEs ()
void	createKDE (const size_t)
void	calculatePIDAKDEMostProbable ()
void	calculatePIDAKDEFullWidthHalfMax ()
void	FillPIDAProperties (unsigned int, unsigned int, unsigned int, anab::Calorimetry const &)

Private Attributes
const float	fPIDA_BOGUS
float	fExponentConstant
float	fMinResRange
float	fMaxResRange
float	fMaxPIDAValue
float	fKDEEvalMaxSigma
float	fKDEEvalStepSize
std::vector< float >	fKDEBandwidths
std::vector< float >	fpida_values
std::vector< float >	fpida_errors
float	fpida_mean
float	fpida_sigma
float	fpida_integral_dedx
float	fpida_integral_pida
std::vector< float >	fpida_kde_mp
std::vector< float >	fpida_kde_fwhm
std::vector< float >	fpida_kde_b
std::vector< std::vector< float > >	fkde_distribution
std::vector< float >	fkde_dist_min
std::vector< float >	fkde_dist_max
util::NormalDistribution	fnormalDist
TTree *	fPIDATree
TH1F *	hPIDAvalues
TH1F *	hPIDAKDE [MAX_BANDWIDTHS]
unsigned int	fPIDAHistNbins
float	fPIDAHistMin
float	fPIDAHistMax
PIDAProperties_t	fPIDAProperties

Detailed Description

Definition at line 48 of file PIDAAlg.h.

Member Typedef Documentation

typedef struct pid::PIDAAlg::PIDAProperties pid::PIDAAlg::PIDAProperties_t

private

Constructor & Destructor Documentation

pid::PIDAAlg::PIDAAlg

(

fhicl::ParameterSet const &

p

)

Title: PIDA Algorithim Class Author: Wes Ketchum (wketc.nosp@m.hum@.nosp@m.lanl..nosp@m.gov), based on ideas/code from Bruce Baller

Description: Algorithm that calculates the PIDA from a calorimetry object Input: anab::Calorimetry Output: PIDA information

Definition at line 24 of file PIDAAlg.cxx.

                                              :
  fPIDA_BOGUS(-9999),
  fExponentConstant(p.get<float>("ExponentConstant",0.42)),
  fMinResRange(p.get<float>("MinResRange",0)),
  fMaxResRange(p.get<float>("MaxResRange",30)),
  fMaxPIDAValue(p.get<float>("MaxPIDAValue",50)),
  fKDEEvalMaxSigma(p.get<float>("KDEEvalMaxSigma",3)),
  fKDEEvalStepSize(p.get<float>("KDEEvalStepSize",0.01)),
  fKDEBandwidths(p.get< std::vector<float> >("KDEBandwidths")),
  fnormalDist(util::NormalDistribution(fKDEEvalMaxSigma,fKDEEvalStepSize)),
  fPIDAHistNbins(p.get<unsigned int>("PIDAHistNbins",100)),
  fPIDAHistMin(p.get<float>("PIDAHistMin",0.0)),
  fPIDAHistMax(p.get<float>("PIDAHistMax",50.0))
{
  ClearInternalData();
}

Member Function Documentation

void pid::PIDAAlg::calculatePIDAIntegral ( std::map< double, double > const &  range_dEdx_map )

private

Definition at line 208 of file PIDAAlg.cxx.

                                                                                 {

  if(range_dEdx_map.size()<2) return;

  fpida_integral_dedx = 0;

  for( std::map<double,double>::const_iterator map_iter = range_dEdx_map.begin();
       map_iter != std::prev(range_dEdx_map.end());
       map_iter++)
    {
      double range_width = std::next(map_iter)->first - map_iter->first;
      fpida_integral_dedx +=  range_width*( std::next(map_iter)->second + 0.5*(map_iter->second-std::next(map_iter)->second));
    }

  fpida_integral_pida = fpida_integral_dedx * (1-fExponentConstant) *
    std::pow( (std::prev(range_dEdx_map.end())->first - range_dEdx_map.begin()->first),(fExponentConstant-1) );
}

void pid::PIDAAlg::calculatePIDAKDEFullWidthHalfMax ( )

private

void pid::PIDAAlg::calculatePIDAKDEMostProbable ( )

private

void pid::PIDAAlg::calculatePIDAMean ( )

private

Definition at line 184 of file PIDAAlg.cxx.

                                  {

  if(fpida_values.size()==0)
    throw "pid::PIDAAlg --- PIDA Values not filled!";

  fpida_mean = 0;
  for(auto const& val : fpida_values)
    fpida_mean += val;
  fpida_mean /= fpida_values.size();

}

void pid::PIDAAlg::calculatePIDASigma ( )

private

Definition at line 196 of file PIDAAlg.cxx.

                                   {

  if(fpida_mean==fPIDA_BOGUS)
    calculatePIDAMean();

  fpida_sigma = 0;
  for(auto const& val : fpida_values)
    fpida_sigma += (fpida_mean-val)*(fpida_mean-val);

  fpida_sigma = std::sqrt(fpida_sigma)/fpida_values.size();
}

void pid::PIDAAlg::ClearInternalData ( )

private

Definition at line 41 of file PIDAAlg.cxx.

                                  {
  fpida_mean = fPIDA_BOGUS;
  fpida_sigma = fPIDA_BOGUS;
  fpida_integral_dedx = fPIDA_BOGUS;
  fpida_integral_pida = fPIDA_BOGUS;
  fpida_values.clear();
  fpida_errors.clear();

  fpida_kde_mp      = std::vector<float>(fKDEBandwidths.size(),fPIDA_BOGUS);
  fpida_kde_fwhm    = std::vector<float>(fKDEBandwidths.size(),fPIDA_BOGUS);
  fpida_kde_b       = std::vector<float>(fKDEBandwidths.size(),fPIDA_BOGUS);
  fkde_distribution = std::vector< std::vector<float> >(fKDEBandwidths.size());
  fkde_dist_min     = std::vector<float>(fKDEBandwidths.size(),fPIDA_BOGUS);
  fkde_dist_max     = std::vector<float>(fKDEBandwidths.size(),fPIDA_BOGUS);
}

void pid::PIDAAlg::createKDE ( const size_t  i_b )

private

Definition at line 226 of file PIDAAlg.cxx.

                                          {

  if(fpida_values.size()==0)
    throw "pid::PIDAAlg --- PIDA Values not filled!";

  //if( fkde_distribution[i_b].size()!=0 ) return;

  if(fKDEBandwidths[i_b]<=0) {
    calculatePIDASigma();
    float bandwidth = fpida_sigma*1.06*std::pow((float)(fpida_values.size()),-0.2);
    fpida_errors = std::vector<float>(fpida_values.size(),bandwidth);
    fpida_kde_b[i_b] = bandwidth;
  }
  else{
    fpida_errors = std::vector<float>(fpida_values.size(),fKDEBandwidths[i_b]);
    fpida_kde_b[i_b] = fKDEBandwidths[i_b];
  }

  const auto min_pida_iterator = std::min_element(fpida_values.begin(),fpida_values.end());
  const size_t min_pida_location = std::distance(fpida_values.begin(),min_pida_iterator);
  fkde_dist_min[i_b] = fpida_values[min_pida_location] - fKDEEvalMaxSigma*fpida_errors[min_pida_location];

  const auto max_pida_iterator = std::max_element(fpida_values.begin(),fpida_values.end());
  const size_t max_pida_location = std::distance(fpida_values.begin(),max_pida_iterator);
  fkde_dist_max[i_b] = fpida_values[max_pida_location] + fKDEEvalMaxSigma*fpida_errors[max_pida_location];

  //make the kde distribution, and get the max value
  const size_t kde_dist_size = (size_t)( (fkde_dist_max[i_b] - fkde_dist_min[i_b])/fKDEEvalStepSize ) + 1;
  fkde_distribution[i_b].resize(kde_dist_size);
  float kde_max=0;
  size_t step_max=0;
  for(size_t i_step=0; i_step<kde_dist_size; i_step++){
    float pida_val = fkde_dist_min[i_b] + i_step*fKDEEvalStepSize;
    fkde_distribution[i_b][i_step]=0;

    for(size_t i_pida=0; i_pida<fpida_values.size(); i_pida++)
      fkde_distribution[i_b][i_step] += fnormalDist.getValue((fpida_values[i_pida]-pida_val)/fpida_errors[i_pida])/fpida_errors[i_pida];

    if(fkde_distribution[i_b][i_step]>kde_max){
      kde_max = fkde_distribution[i_b][i_step];
      step_max = i_step;
      fpida_kde_mp[i_b] = pida_val;
    }
  }

  //now get fwhm
  float half_max = 0.5*fpida_kde_mp[i_b];
  float low_width=0;
  for(size_t i_step=step_max; i_step>0; i_step--){
    if(fkde_distribution[i_b][i_step] < half_max) break;
    low_width += fKDEEvalStepSize;
  }
  float high_width=0;
  for(size_t i_step=step_max; i_step<kde_dist_size; i_step++){
    if(fkde_distribution[i_b][i_step] < half_max) break;
    high_width += fKDEEvalStepSize;
  }
  fpida_kde_fwhm[i_b] = low_width+high_width;

}

void pid::PIDAAlg::createKDEs ( )

private

Definition at line 287 of file PIDAAlg.cxx.

                           {
  for(size_t i_b=0; i_b < fKDEBandwidths.size(); i_b++)
    createKDE(i_b);
}

void pid::PIDAAlg::FillPIDAProperties ( unsigned int  run, unsigned int  event, unsigned int  calo_index, anab::Calorimetry const &  calo  )

private

Definition at line 292 of file PIDAAlg.cxx.

                                                                  {

  fPIDAProperties.run = run;
  fPIDAProperties.event = event;
  fPIDAProperties.calo_index = calo_index;
  fPIDAProperties.planeid = calo.PlaneID().Plane;
  fPIDAProperties.trk_range = calo.Range();
  fPIDAProperties.calo_KE = calo.KineticEnergy();

  fPIDAProperties.n_bandwidths = fKDEBandwidths.size();
  for(size_t i_b=0; i_b<fPIDAProperties.n_bandwidths; i_b++)

  calculatePIDASigma();
  fPIDAProperties.n_pid_pts = fpida_values.size();
  fPIDAProperties.mean = fpida_mean;
  fPIDAProperties.sigma = fpida_sigma;

  fPIDAProperties.integral_dedx = fpida_integral_dedx;
  fPIDAProperties.integral_pida = fpida_integral_pida;

  createKDEs();
  for(size_t i_b=0; i_b<fPIDAProperties.n_bandwidths; i_b++){
    fPIDAProperties.kde_mp[i_b]        = fpida_kde_mp[i_b];
    fPIDAProperties.kde_fwhm[i_b]      = fpida_kde_fwhm[i_b];
    fPIDAProperties.kde_bandwidth[i_b] = fpida_kde_b[i_b];
  }

  hPIDAvalues->Reset();
  for(auto const& val: fpida_values)
    hPIDAvalues->Fill(val);

  for(size_t i_b=0; i_b<fPIDAProperties.n_bandwidths; i_b++){
    hPIDAKDE[i_b]->Reset();
    for(size_t i_step=0; i_step<fkde_distribution[i_b].size(); i_step++)
      hPIDAKDE[i_b]->AddBinContent(hPIDAKDE[i_b]->FindBin(i_step*fKDEEvalStepSize+fkde_dist_min[i_b]),
                                   fkde_distribution[i_b][i_step]);
  }

  fPIDATree->Fill();
}

void pid::PIDAAlg::FillPIDATree	(	unsigned int	run,
		unsigned int	event,
		unsigned int	calo_index,
		anab::Calorimetry const &	calo
	)

Definition at line 176 of file PIDAAlg.cxx.

                                                            {
  RunPIDAAlg(calo);
  FillPIDAProperties(run,event,calo_index,calo);
}

float pid::PIDAAlg::getKDEBandwidth ( const size_t  i_b )

inline

Definition at line 63 of file PIDAAlg.h.

{ return fKDEBandwidths.at(i_b); }

size_t pid::PIDAAlg::getNKDEBandwidths ( )

inline

Definition at line 62 of file PIDAAlg.h.

{ return fKDEBandwidths.size(); }

const std::vector< float > & pid::PIDAAlg::getPIDAErrors

(

)

Definition at line 127 of file PIDAAlg.cxx.

                                                 {
  return fpida_errors;
}

float pid::PIDAAlg::getPIDAKDEFullWidthHalfMax

(

const size_t

i_b

)

Definition at line 116 of file PIDAAlg.cxx.

                                                            {
  if(fpida_kde_fwhm[i_b]==fPIDA_BOGUS)
    createKDE(i_b);

  return fpida_kde_fwhm[i_b];
}

float pid::PIDAAlg::getPIDAKDEMostProbable

(

const size_t

i_b

)

Definition at line 109 of file PIDAAlg.cxx.

                                                        {
  if(fpida_kde_mp[i_b]==fPIDA_BOGUS)
    createKDE(i_b);

  return fpida_kde_mp[i_b];
}

float pid::PIDAAlg::getPIDAMean

(

)

Definition at line 95 of file PIDAAlg.cxx.

                             {
  if(fpida_mean==fPIDA_BOGUS)
    calculatePIDAMean();

  return fpida_mean;
}

float pid::PIDAAlg::getPIDASigma

(

)

Definition at line 102 of file PIDAAlg.cxx.

                              {
  if(fpida_sigma==fPIDA_BOGUS)
    calculatePIDASigma();

  return fpida_sigma;
}

const std::vector< float > & pid::PIDAAlg::getPIDAValues

(

)

Definition at line 123 of file PIDAAlg.cxx.

                                                 {
  return fpida_values;
}

void pid::PIDAAlg::PrintPIDAValues

(

)

Definition at line 336 of file PIDAAlg.cxx.

                                {
  for(size_t i_pida=0; i_pida<fpida_values.size(); i_pida++)
    std::cout << "\tPIDA --- " << i_pida << "\t" << fpida_values[i_pida] << std::endl;
}

void pid::PIDAAlg::RunPIDAAlg	(	std::vector< float > const &	resRange,
		std::vector< float > const &	dEdx
	)

Definition at line 146 of file PIDAAlg.cxx.

                                                           {

  ClearInternalData();

  fpida_values.reserve( resRange.size() );
  fpida_errors.reserve( resRange.size() );

  std::map<double,double> range_dEdx_map;

  for(size_t i_r=0; i_r<resRange.size(); i_r++){
    if(resRange[i_r]>fMaxResRange || resRange[i_r]<fMinResRange) continue;

    range_dEdx_map[ resRange[i_r] ] = dEdx[i_r];

    float val = dEdx[i_r]*std::pow(resRange[i_r],fExponentConstant);
    if(val < fMaxPIDAValue){
      fpida_values.push_back(val);
      //fpida_errors.push_back(0);
    }

  }

  calculatePIDAIntegral(range_dEdx_map);

  if(fpida_values.size()==0)
    fpida_values.push_back(-99);

}

void pid::PIDAAlg::RunPIDAAlg

(

anab::Calorimetry const &

calo

)

Definition at line 131 of file PIDAAlg.cxx.

                                                      {
  std::vector<float> const& resRange = calo.ResidualRange();
  std::vector<float> const& dEdx     = calo.dEdx();
  RunPIDAAlg(resRange,dEdx);
}

void pid::PIDAAlg::RunPIDAAlg	(	anab::Calorimetry const &	calo,
		float &	mean,
		float &	sigma
	)

Definition at line 137 of file PIDAAlg.cxx.

{
  RunPIDAAlg(calo);
  mean = getPIDAMean();
  sigma = getPIDASigma();
}

void pid::PIDAAlg::setExponentConstant ( float const &  ex )

inline

Definition at line 69 of file PIDAAlg.h.

{ fExponentConstant = ex; }

void pid::PIDAAlg::SetPIDATree	(	TTree *	tree,
		TH1F *	hist_vals,
		std::vector< TH1F * >	hist_kde
	)

Definition at line 57 of file PIDAAlg.cxx.

                                                                                   {

  if(hist_kde.size()>MAX_BANDWIDTHS)
    throw "Error: input histograms larger than max allowed bandwidths.";
  if(hist_kde.size()!=fKDEBandwidths.size())
    throw "Error: input histograms do not have same size as bandwidths.";

  fPIDATree = tree;

  hPIDAvalues = hist_vals;
  hPIDAvalues->SetNameTitle("hPIDAvalues","PIDA Distribution");
  hPIDAvalues->SetBins(fPIDAHistNbins,fPIDAHistMin,fPIDAHistMax);

  for(size_t i_hist=0; i_hist<hist_kde.size(); i_hist++){
    hPIDAKDE[i_hist] = hist_kde[i_hist];

    std::stringstream hname,htitle;
    hname << "hPIDAKDE_" << i_hist;
    htitle << "PIDA KDE-smoothed Distribution, Bandwidth=" << fKDEBandwidths.at(i_hist);

    hPIDAKDE[i_hist]->SetNameTitle(hname.str().c_str(),htitle.str().c_str());
    hPIDAKDE[i_hist]->SetBins(fPIDAHistNbins,fPIDAHistMin,fPIDAHistMax);
  }

  fPIDATree->Branch("pida",&fPIDAProperties,fPIDAProperties.leaf_structure.c_str());
  fPIDATree->Branch("hpida_vals","TH1F",hPIDAvalues);
  fPIDATree->Branch("n_bandwidths",&(fPIDAProperties.n_bandwidths),"n_bandwidths/i");
  fPIDATree->Branch("kde_bandwidth",fPIDAProperties.kde_bandwidth,"kde_bandwidth[n_bandwidths]/F");
  fPIDATree->Branch("kde_mp",fPIDAProperties.kde_mp,"kde_mp[n_bandwidths]/F");
  fPIDATree->Branch("kde_fwhm",fPIDAProperties.kde_fwhm,"kde_fwhm[n_bandwidths]/F");
  for(size_t i_hist=0; i_hist<hist_kde.size(); i_hist++){
    std::stringstream bname;
    bname << "hpida_kde_" << i_hist;
    fPIDATree->Branch(bname.str().c_str(),"TH1F",hPIDAKDE[i_hist]);
  }

}

Member Data Documentation

float pid::PIDAAlg::fExponentConstant

private

Definition at line 78 of file PIDAAlg.h.

std::vector<float> pid::PIDAAlg::fkde_dist_max

private

Definition at line 110 of file PIDAAlg.h.

std::vector<float> pid::PIDAAlg::fkde_dist_min

private

Definition at line 109 of file PIDAAlg.h.

std::vector< std::vector<float> > pid::PIDAAlg::fkde_distribution

private

Definition at line 108 of file PIDAAlg.h.

std::vector<float> pid::PIDAAlg::fKDEBandwidths

private

Definition at line 84 of file PIDAAlg.h.

float pid::PIDAAlg::fKDEEvalMaxSigma

private

Definition at line 82 of file PIDAAlg.h.

float pid::PIDAAlg::fKDEEvalStepSize

private

Definition at line 83 of file PIDAAlg.h.

float pid::PIDAAlg::fMaxPIDAValue

private

Definition at line 81 of file PIDAAlg.h.

float pid::PIDAAlg::fMaxResRange

private

Definition at line 80 of file PIDAAlg.h.

float pid::PIDAAlg::fMinResRange

private

Definition at line 79 of file PIDAAlg.h.

util::NormalDistribution pid::PIDAAlg::fnormalDist

private

Definition at line 112 of file PIDAAlg.h.

const float pid::PIDAAlg::fPIDA_BOGUS

private

Definition at line 76 of file PIDAAlg.h.

std::vector<float> pid::PIDAAlg::fpida_errors

private

Definition at line 87 of file PIDAAlg.h.

float pid::PIDAAlg::fpida_integral_dedx

private

Definition at line 90 of file PIDAAlg.h.

float pid::PIDAAlg::fpida_integral_pida

private

Definition at line 91 of file PIDAAlg.h.

std::vector<float> pid::PIDAAlg::fpida_kde_b

private

Definition at line 105 of file PIDAAlg.h.

std::vector<float> pid::PIDAAlg::fpida_kde_fwhm

private

Definition at line 104 of file PIDAAlg.h.

std::vector<float> pid::PIDAAlg::fpida_kde_mp

private

Definition at line 103 of file PIDAAlg.h.

float pid::PIDAAlg::fpida_mean

private

Definition at line 88 of file PIDAAlg.h.

float pid::PIDAAlg::fpida_sigma

private

Definition at line 89 of file PIDAAlg.h.

std::vector<float> pid::PIDAAlg::fpida_values

private

Definition at line 86 of file PIDAAlg.h.

float pid::PIDAAlg::fPIDAHistMax

private

Definition at line 119 of file PIDAAlg.h.

float pid::PIDAAlg::fPIDAHistMin

private

Definition at line 118 of file PIDAAlg.h.

unsigned int pid::PIDAAlg::fPIDAHistNbins

private

Definition at line 117 of file PIDAAlg.h.

PIDAProperties_t pid::PIDAAlg::fPIDAProperties

private

Definition at line 144 of file PIDAAlg.h.

TTree* pid::PIDAAlg::fPIDATree

private

Definition at line 114 of file PIDAAlg.h.

TH1F* pid::PIDAAlg::hPIDAKDE[MAX_BANDWIDTHS]

private

Definition at line 116 of file PIDAAlg.h.

TH1F* pid::PIDAAlg::hPIDAvalues

private

Definition at line 115 of file PIDAAlg.h.

---

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

• larana/larana/ParticleIdentification/PIDAAlg.h
• larana/larana/ParticleIdentification/PIDAAlg.cxx