pid::Chi2PIDAlg Class Reference

#include "Chi2PIDAlg.h"

Public Member Functions
	Chi2PIDAlg (fhicl::ParameterSet const &pset)
std::bitset< 8 >	GetBitset (geo::PlaneID planeID)
anab::ParticleID	DoParticleID (const std::vector< art::Ptr< anab::Calorimetry >> &calo)

Private Attributes
std::string	fTemplateFile
bool	fUseMedian
std::string	fROOTfile
TProfile *	dedx_range_pro
	proton template More...
TProfile *	dedx_range_ka
	kaon template More...
TProfile *	dedx_range_pi
	pion template More...
TProfile *	dedx_range_mu
	muon template More...

Detailed Description

Definition at line 30 of file Chi2PIDAlg.h.

Constructor & Destructor Documentation

pid::Chi2PIDAlg::Chi2PIDAlg

(

fhicl::ParameterSet const &

pset

)

Definition at line 28 of file Chi2PIDAlg.cxx.

References dedx_range_ka, dedx_range_mu, dedx_range_pi, dedx_range_pro, file, fROOTfile, fTemplateFile, fUseMedian, and fhicl::ParameterSet::get().

{
  fTemplateFile = pset.get<std::string>("TemplateFile");
  fUseMedian = pset.get<bool>("UseMedian");
  //fCalorimetryModuleLabel = pset.get< std::string >("CalorimetryModuleLabel");

  cet::search_path sp("FW_SEARCH_PATH");

  if (!sp.find_file(fTemplateFile, fROOTfile))
    throw cet::exception("Chi2ParticleID") << "cannot find the root template file: \n"
                                           << fTemplateFile << "\n bail ungracefully.\n";
  TFile* file = TFile::Open(fROOTfile.c_str());
  dedx_range_pro = (TProfile*)file->Get("dedx_range_pro");
  dedx_range_ka = (TProfile*)file->Get("dedx_range_ka");
  dedx_range_pi = (TProfile*)file->Get("dedx_range_pi");
  dedx_range_mu = (TProfile*)file->Get("dedx_range_mu");

  //  std::cout<<"Chi2PIDAlg configuration:"<<std::endl;
  //  std::cout<<"Template file: "<<fROOTfile<<std::endl;
  //  std::cout<<"fUseMedian: "<<fUseMedian<<std::endl;
}

Member Function Documentation

anab::ParticleID pid::Chi2PIDAlg::DoParticleID

(

const std::vector< art::Ptr< anab::Calorimetry >> &

calo

)

Definition at line 62 of file Chi2PIDAlg.cxx.

References bin, anab::Calorimetry::DeadWireResRC(), anab::Calorimetry::dEdx(), dedx_range_ka, dedx_range_mu, dedx_range_pi, dedx_range_pro, e, anab::sParticleIDAlgScores::fAlgName, anab::sParticleIDAlgScores::fAssumedPdg, anab::sParticleIDAlgScores::fNdf, anab::sParticleIDAlgScores::fPlaneMask, anab::sParticleIDAlgScores::fTrackDir, fUseMedian, anab::sParticleIDAlgScores::fValue, anab::sParticleIDAlgScores::fVariableType, GetBitset(), anab::kForward, anab::kGOF, anab::kPIDA, anab::Calorimetry::PlaneID(), and anab::Calorimetry::ResidualRange().

Referenced by pid::Chi2ParticleID::produce().

{

  std::vector<anab::sParticleIDAlgScores> AlgScoresVec;
  geo::PlaneID plid;

  for (size_t i_calo = 0; i_calo < calos.size(); i_calo++) {

    art::Ptr<anab::Calorimetry> calo = calos.at(i_calo);
    if (i_calo == 0)
      plid = calo->PlaneID();
    else if (plid != calo->PlaneID())
      throw cet::exception("Chi2PIDAlg") << "PlaneID mismatch: " << plid << ", " << calo->PlaneID();
    int npt = 0;
    double chi2pro = 0;
    double chi2ka = 0;
    double chi2pi = 0;
    double chi2mu = 0;
    double PIDA = 0; //by Bruce Baller
    std::vector<double> vpida;
    std::vector<float> trkdedx = calo->dEdx();
    std::vector<float> trkres = calo->ResidualRange();
    std::vector<float> deadwireresrc = calo->DeadWireResRC();

    int used_trkres = 0;
    int nbins_dedx_range = dedx_range_pro->GetNbinsX();
    for (unsigned i = 0; i < trkdedx.size(); ++i) { //hits
      //ignore the first and the last point
      if (i == 0 || i == trkdedx.size() - 1) continue;
      if (trkres[i] < 30) {
        double PIDAi = trkdedx[i] * pow(trkres[i], 0.42);
        PIDA += PIDAi;
        vpida.push_back(PIDAi);
        used_trkres++;
      }
      if (trkdedx[i] > 1000) continue; //protect against large pulse height
      int bin = dedx_range_pro->FindBin(trkres[i]);
      if (bin >= 1 && bin <= nbins_dedx_range) {
        double bincpro = dedx_range_pro->GetBinContent(bin);
        if (bincpro < 1e-6) { //for 0 bin content, using neighboring bins
          bincpro =
            (dedx_range_pro->GetBinContent(bin - 1) + dedx_range_pro->GetBinContent(bin + 1)) / 2;
        }
        double bincka = dedx_range_ka->GetBinContent(bin);
        if (bincka < 1e-6) {
          bincka =
            (dedx_range_ka->GetBinContent(bin - 1) + dedx_range_ka->GetBinContent(bin + 1)) / 2;
        }
        double bincpi = dedx_range_pi->GetBinContent(bin);
        if (bincpi < 1e-6) {
          bincpi =
            (dedx_range_pi->GetBinContent(bin - 1) + dedx_range_pi->GetBinContent(bin + 1)) / 2;
        }
        double bincmu = dedx_range_mu->GetBinContent(bin);
        if (bincmu < 1e-6) {
          bincmu =
            (dedx_range_mu->GetBinContent(bin - 1) + dedx_range_mu->GetBinContent(bin + 1)) / 2;
        }
        double binepro = dedx_range_pro->GetBinError(bin);
        if (binepro < 1e-6) {
          binepro =
            (dedx_range_pro->GetBinError(bin - 1) + dedx_range_pro->GetBinError(bin + 1)) / 2;
        }
        double bineka = dedx_range_ka->GetBinError(bin);
        if (bineka < 1e-6) {
          bineka = (dedx_range_ka->GetBinError(bin - 1) + dedx_range_ka->GetBinError(bin + 1)) / 2;
        }
        double binepi = dedx_range_pi->GetBinError(bin);
        if (binepi < 1e-6) {
          binepi = (dedx_range_pi->GetBinError(bin - 1) + dedx_range_pi->GetBinError(bin + 1)) / 2;
        }
        double binemu = dedx_range_mu->GetBinError(bin);
        if (binemu < 1e-6) {
          binemu = (dedx_range_mu->GetBinError(bin - 1) + dedx_range_mu->GetBinError(bin + 1)) / 2;
        }
        //double errke = 0.05*trkdedx[i];   //5% KE resolution
        double errdedx = 0.04231 + 0.0001783 * trkdedx[i] * trkdedx[i]; //resolution on dE/dx
        errdedx *= trkdedx[i];

        double errdedx_square = errdedx * errdedx;
        chi2pro += cet::square(trkdedx[i] - bincpro) / (binepro * binepro + errdedx_square);
        chi2ka += cet::square(trkdedx[i] - bincka) / (bineka * bineka + errdedx_square);
        chi2pi += cet::square(trkdedx[i] - bincpi) / (binepi * binepi + errdedx_square);
        chi2mu += cet::square(trkdedx[i] - bincmu) / (binemu * binemu + errdedx_square);

        //std::cout<<i<<" "<<trkdedx[i]<<" "<<trkres[i]<<" "<<bincpro<<std::endl;
        ++npt;
      }
    }

    anab::sParticleIDAlgScores chi2proton;
    anab::sParticleIDAlgScores chi2kaon;
    anab::sParticleIDAlgScores chi2pion;
    anab::sParticleIDAlgScores chi2muon;
    anab::sParticleIDAlgScores pida_mean;
    anab::sParticleIDAlgScores pida_median;

    //anab::ParticleID pidOut;
    if (npt) {

      chi2proton.fAlgName = "Chi2";
      chi2proton.fVariableType = anab::kGOF;
      chi2proton.fTrackDir = anab::kForward;
      chi2proton.fAssumedPdg = 2212;
      chi2proton.fPlaneMask = GetBitset(calo->PlaneID());
      chi2proton.fNdf = npt;
      chi2proton.fValue = chi2pro / npt;

      chi2muon.fAlgName = "Chi2";
      chi2muon.fVariableType = anab::kGOF;
      chi2muon.fTrackDir = anab::kForward;
      chi2muon.fAssumedPdg = 13;
      chi2muon.fPlaneMask = GetBitset(calo->PlaneID());
      chi2muon.fNdf = npt;
      chi2muon.fValue = chi2mu / npt;

      chi2kaon.fAlgName = "Chi2";
      chi2kaon.fVariableType = anab::kGOF;
      chi2kaon.fTrackDir = anab::kForward;
      chi2kaon.fAssumedPdg = 321;
      chi2kaon.fPlaneMask = GetBitset(calo->PlaneID());
      chi2kaon.fNdf = npt;
      chi2kaon.fValue = chi2ka / npt;

      chi2pion.fAlgName = "Chi2";
      chi2pion.fVariableType = anab::kGOF;
      chi2pion.fTrackDir = anab::kForward;
      chi2pion.fAssumedPdg = 211;
      chi2pion.fPlaneMask = GetBitset(calo->PlaneID());
      chi2pion.fNdf = npt;
      chi2pion.fValue = chi2pi / npt;

      AlgScoresVec.push_back(chi2proton);
      AlgScoresVec.push_back(chi2muon);
      AlgScoresVec.push_back(chi2kaon);
      AlgScoresVec.push_back(chi2pion);
    }

    //if (trkdedx.size()) pidOut.fPIDA = PIDA/trkdedx.size();
    if (used_trkres > 0) {
      if (fUseMedian) {
        pida_median.fAlgName = "PIDA_median";
        pida_median.fVariableType = anab::kPIDA;
        pida_median.fTrackDir = anab::kForward;
        pida_median.fValue = TMath::Median(vpida.size(), &vpida[0]);
        pida_median.fPlaneMask = GetBitset(calo->PlaneID());
        AlgScoresVec.push_back(pida_median);
      }
      else { // use mean
        pida_mean.fAlgName = "PIDA_mean";
        pida_mean.fVariableType = anab::kPIDA;
        pida_mean.fTrackDir = anab::kForward;
        pida_mean.fValue = PIDA / used_trkres;
        pida_mean.fPlaneMask = GetBitset(calo->PlaneID());
        AlgScoresVec.push_back(pida_mean);
      }
    }
  }

  anab::ParticleID pidOut(AlgScoresVec, plid);

  return pidOut;
}

std::bitset< 8 > pid::Chi2PIDAlg::GetBitset

(

geo::PlaneID

planeID

)

Helper function to go from geo::PlaneID to a bitset

Definition at line 51 of file Chi2PIDAlg.cxx.

References geo::PlaneID::Plane.

Referenced by DoParticleID().

{

  std::bitset<8> thisBitset;

  thisBitset.set(planeID.Plane);

  return thisBitset;
}

Member Data Documentation

TProfile* pid::Chi2PIDAlg::dedx_range_ka

private

kaon template

Definition at line 49 of file Chi2PIDAlg.h.

Referenced by Chi2PIDAlg(), and DoParticleID().

TProfile* pid::Chi2PIDAlg::dedx_range_mu

private

muon template

Definition at line 51 of file Chi2PIDAlg.h.

Referenced by Chi2PIDAlg(), and DoParticleID().

TProfile* pid::Chi2PIDAlg::dedx_range_pi

private

pion template

Definition at line 50 of file Chi2PIDAlg.h.

Referenced by Chi2PIDAlg(), and DoParticleID().

TProfile* pid::Chi2PIDAlg::dedx_range_pro

private

proton template

Definition at line 48 of file Chi2PIDAlg.h.

Referenced by Chi2PIDAlg(), and DoParticleID().

std::string pid::Chi2PIDAlg::fROOTfile

private

Definition at line 46 of file Chi2PIDAlg.h.

Referenced by Chi2PIDAlg().

std::string pid::Chi2PIDAlg::fTemplateFile

private

Definition at line 43 of file Chi2PIDAlg.h.

Referenced by Chi2PIDAlg().

bool pid::Chi2PIDAlg::fUseMedian

private

Definition at line 44 of file Chi2PIDAlg.h.

Referenced by Chi2PIDAlg(), and DoParticleID().

---

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

• larana/v09_15_05/source/larana/ParticleIdentification/Chi2PIDAlg.h
• larana/v09_15_05/source/larana/ParticleIdentification/Chi2PIDAlg.cxx