LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
Public Member Functions | Private Attributes | List of all members
pid::Chi2PIDAlg Class Reference

#include "Chi2PIDAlg.h"

Public Member Functions

 Chi2PIDAlg (fhicl::ParameterSet const &pset)
 
virtual ~Chi2PIDAlg ()
 
void reconfigure (fhicl::ParameterSet const &pset)
 
void DoParticleID (art::Ptr< anab::Calorimetry > calo, anab::ParticleID &pidOut)
 

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 28 of file Chi2PIDAlg.h.

Constructor & Destructor Documentation

pid::Chi2PIDAlg::Chi2PIDAlg ( fhicl::ParameterSet const &  pset)

Definition at line 35 of file Chi2PIDAlg.cxx.

References reconfigure().

36 {
37  this->reconfigure(pset);
38 }
pid::Chi2PIDAlg::reconfigure
void reconfigure(fhicl::ParameterSet const &pset)
Definition: Chi2PIDAlg.cxx:46
pid::Chi2PIDAlg::~Chi2PIDAlg ( )
virtual

Definition at line 41 of file Chi2PIDAlg.cxx.

42 {
43 }

Member Function Documentation

void pid::Chi2PIDAlg::DoParticleID ( art::Ptr< anab::Calorimetry >  calo,
anab::ParticleID &  pidOut 
)

Definition at line 73 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::ParticleID::fChi2Kaon, anab::ParticleID::fChi2Muon, anab::ParticleID::fChi2Pion, anab::ParticleID::fChi2Proton, anab::ParticleID::fDeltaChi2, anab::ParticleID::fMinChi2, anab::ParticleID::fMissingE, anab::ParticleID::fMissingEavg, anab::ParticleID::fNdf, anab::ParticleID::fPdg, anab::ParticleID::fPIDA, anab::ParticleID::fPlaneID, fUseMedian, anab::Calorimetry::PlaneID(), anab::Calorimetry::ResidualRange(), and anab::Calorimetry::TrkPitchC().

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

74  {
75  int npt = 0;
76  double chi2pro = 0;
77  double chi2ka = 0;
78  double chi2pi = 0;
79  double chi2mu = 0;
80  double trkpitchc = calo->TrkPitchC();
81  double avgdedx = 0;
82  double PIDA = 0; //by Bruce Baller
83  std::vector<double> vpida;
84  std::vector<float> trkdedx = calo->dEdx();
85  std::vector<float> trkres = calo->ResidualRange();
86  std::vector<float> deadwireresrc = calo->DeadWireResRC();
87  pidOut.fPlaneID = calo->PlaneID();
88 
89  int used_trkres = 0;
90  for (unsigned i = 0; i<trkdedx.size(); ++i){//hits
91  //ignore the first and the last point
92  if (i==0 || i==trkdedx.size()-1) continue;
93  avgdedx += trkdedx[i];
94  if(trkres[i] < 30) {
95  PIDA += trkdedx[i]*pow(trkres[i],0.42);
96  vpida.push_back(trkdedx[i]*pow(trkres[i],0.42));
97  used_trkres++;
98  }
99  if (trkdedx[i]>1000) continue; //protect against large pulse height
100  int bin = dedx_range_pro->FindBin(trkres[i]);
101  if (bin>=1&&bin<=dedx_range_pro->GetNbinsX()){
102  double bincpro = dedx_range_pro->GetBinContent(bin);
103  if (bincpro<1e-6){//for 0 bin content, using neighboring bins
104  bincpro = (dedx_range_pro->GetBinContent(bin-1)+dedx_range_pro->GetBinContent(bin+1))/2;
105  }
106  double bincka = dedx_range_ka->GetBinContent(bin);
107  if (bincka<1e-6){
108  bincka = (dedx_range_ka->GetBinContent(bin-1)+dedx_range_ka->GetBinContent(bin+1))/2;
109  }
110  double bincpi = dedx_range_pi->GetBinContent(bin);
111  if (bincpi<1e-6){
112  bincpi = (dedx_range_pi->GetBinContent(bin-1)+dedx_range_pi->GetBinContent(bin+1))/2;
113  }
114  double bincmu = dedx_range_mu->GetBinContent(bin);
115  if (bincmu<1e-6){
116  bincmu = (dedx_range_mu->GetBinContent(bin-1)+dedx_range_mu->GetBinContent(bin+1))/2;
117  }
118  double binepro = dedx_range_pro->GetBinError(bin);
119  if (binepro<1e-6){
120  binepro = (dedx_range_pro->GetBinError(bin-1)+dedx_range_pro->GetBinError(bin+1))/2;
121  }
122  double bineka = dedx_range_ka->GetBinError(bin);
123  if (bineka<1e-6){
124  bineka = (dedx_range_ka->GetBinError(bin-1)+dedx_range_ka->GetBinError(bin+1))/2;
125  }
126  double binepi = dedx_range_pi->GetBinError(bin);
127  if (binepi<1e-6){
128  binepi = (dedx_range_pi->GetBinError(bin-1)+dedx_range_pi->GetBinError(bin+1))/2;
129  }
130  double binemu = dedx_range_mu->GetBinError(bin);
131  if (binemu<1e-6){
132  binemu = (dedx_range_mu->GetBinError(bin-1)+dedx_range_mu->GetBinError(bin+1))/2;
133  }
134  //double errke = 0.05*trkdedx[i]; //5% KE resolution
135  double errdedx = 0.04231+0.0001783*trkdedx[i]*trkdedx[i]; //resolution on dE/dx
136  errdedx *= trkdedx[i];
137  chi2pro += pow((trkdedx[i]-bincpro)/std::sqrt(pow(binepro,2)+pow(errdedx,2)),2);
138  chi2ka += pow((trkdedx[i]-bincka)/std::sqrt(pow(bineka,2)+pow(errdedx,2)),2);
139  chi2pi += pow((trkdedx[i]-bincpi)/std::sqrt(pow(binepi,2)+pow(errdedx,2)),2);
140  chi2mu += pow((trkdedx[i]-bincmu)/std::sqrt(pow(binemu,2)+pow(errdedx,2)),2);
141  //std::cout<<i<<" "<<trkdedx[i]<<" "<<trkres[i]<<" "<<bincpro<<std::endl;
142  ++npt;
143  }
144  }
145 
146  //anab::ParticleID pidOut;
147  if (npt){
148  pidOut.fNdf = npt;
149  pidOut.fChi2Proton = chi2pro/npt;
150  pidOut.fChi2Kaon = chi2ka/npt;
151  pidOut.fChi2Pion = chi2pi/npt;
152  pidOut.fChi2Muon = chi2mu/npt;
153  double chi2[4] = {chi2pro/npt,chi2ka/npt,chi2pi/npt,chi2mu/npt};
154  double pdg[4] = {2212,321,211,13};
155  double chi2min = 1e20;
156  int imin = -1;
157  double chi2min2 = 1e20;
158  //int imin2;
159  // find the minimal chi2 and next-to-minimal chi2
160  for (int ichi2 = 0; ichi2<4; ++ichi2){
161  if (chi2[ichi2]<chi2min){
162  imin = ichi2;
163  chi2min2 = chi2min;
164  chi2min = chi2[ichi2];
165  }
166  else if (chi2[ichi2]<chi2min2){
167  //imin2 = ichi2;
168  chi2min2 = chi2[ichi2];
169  }
170  }
171  if (imin>-1){
172  pidOut.fPdg = pdg[imin];
173  pidOut.fMinChi2 = chi2min;
174  pidOut.fDeltaChi2 = chi2min2 - chi2min;
175  }
176  }
177  //if (trkdedx.size()) pidOut.fPIDA = PIDA/trkdedx.size();
178  if(used_trkres > 0){
179  if (fUseMedian){
180  pidOut.fPIDA = TMath::Median(vpida.size(), &vpida[0]);
181  }
182  else{//use mean
183  pidOut.fPIDA = PIDA/used_trkres;
184  }
185  }
186  double missinge = 0;
187  double missingeavg = 0;
188  for (unsigned i = 0; i<deadwireresrc.size(); ++i){
189  int bin = dedx_range_pro->FindBin(deadwireresrc[i]);
190  //std::cout<<i<<" "<<deadwireresrc[i]<<" "<<bin<<std::endl;
191  if (bin<1) continue;
192  if (bin>dedx_range_pro->GetNbinsX()) bin = dedx_range_pro->GetNbinsX();
193  if (pidOut.fPdg==2212){
194  missinge += dedx_range_pro->GetBinContent(bin)*trkpitchc;
195  }
196  else if (pidOut.fPdg==321){
197  missinge += dedx_range_ka->GetBinContent(bin)*trkpitchc;
198  }
199  else if (pidOut.fPdg==211){
200  missinge += dedx_range_pi->GetBinContent(bin)*trkpitchc;
201  }
202  else if (pidOut.fPdg==13){
203  missinge += dedx_range_mu->GetBinContent(bin)*trkpitchc;
204  }
205  //std::cout<<bin<<" "<<dedx_range_pro->GetBinContent(bin)*trkpitchc<<std::endl;
206  }
207  if (trkdedx.size()) missingeavg = avgdedx/trkdedx.size()*trkpitchc*deadwireresrc.size();
208  //std::cout<<trkIter<<" "<<pid<<std::endl;
209  pidOut.fMissingE = missinge;
210  pidOut.fMissingEavg = missingeavg;
211 }
pid::Chi2PIDAlg::dedx_range_mu
TProfile * dedx_range_mu
muon template
Definition: Chi2PIDAlg.h:50
anab::ParticleID::fMissingEavg
double fMissingEavg
missing energy from dead wires using average dEdx
Definition: ParticleID.h:34
anab::ParticleID::fMinChi2
double fMinChi2
Minimum reduced chi2.
Definition: ParticleID.h:27
anab::Calorimetry::PlaneID
const geo::PlaneID & PlaneID() const
Definition: Calorimetry.h:119
pid::Chi2PIDAlg::dedx_range_pro
TProfile * dedx_range_pro
proton template
Definition: Chi2PIDAlg.h:47
pid::Chi2PIDAlg::dedx_range_pi
TProfile * dedx_range_pi
pion template
Definition: Chi2PIDAlg.h:49
anab::ParticleID::fChi2Muon
double fChi2Muon
reduced chi2 using muon template
Definition: ParticleID.h:32
anab::ParticleID::fNdf
int fNdf
ndf for chi2 test
Definition: ParticleID.h:26
anab::Calorimetry::DeadWireResRC
const std::vector< float > & DeadWireResRC() const
Definition: Calorimetry.h:107
anab::Calorimetry::ResidualRange
const std::vector< float > & ResidualRange() const
Definition: Calorimetry.h:106
anab::ParticleID::fDeltaChi2
double fDeltaChi2
difference between two lowest reduced chi2&#39;s
Definition: ParticleID.h:28
bin
float bin[41]
Definition: plottest35.C:14
anab::Calorimetry::dEdx
const std::vector< float > & dEdx() const
Definition: Calorimetry.h:104
anab::Calorimetry::TrkPitchC
float TrkPitchC() const
Definition: Calorimetry.h:111
anab::ParticleID::fChi2Kaon
double fChi2Kaon
reduced chi2 using kaon template
Definition: ParticleID.h:30
anab::ParticleID::fPdg
int fPdg
determined particle ID
Definition: ParticleID.h:25
anab::ParticleID::fChi2Proton
double fChi2Proton
reduced chi2 using proton template
Definition: ParticleID.h:29
anab::ParticleID::fPIDA
double fPIDA
PID developed by Bruce Baller.
Definition: ParticleID.h:35
anab::ParticleID::fPlaneID
geo::PlaneID fPlaneID
Definition: ParticleID.h:36
pid::Chi2PIDAlg::dedx_range_ka
TProfile * dedx_range_ka
kaon template
Definition: Chi2PIDAlg.h:48
e
Float_t e
Definition: plot.C:34
anab::ParticleID::fChi2Pion
double fChi2Pion
reduced chi2 using pion template
Definition: ParticleID.h:31
anab::ParticleID::fMissingE
double fMissingE
missing energy from dead wires for contained particle
Definition: ParticleID.h:33
void pid::Chi2PIDAlg::reconfigure ( fhicl::ParameterSet const &  pset)

Definition at line 46 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().

Referenced by Chi2PIDAlg(), and pid::Chi2ParticleID::reconfigure().

47 {
48  fTemplateFile = pset.get< std::string >("TemplateFile");
49  fUseMedian = pset.get< bool >("UseMedian");
50  //fCalorimetryModuleLabel = pset.get< std::string >("CalorimetryModuleLabel");
51 
52  cet::search_path sp("FW_SEARCH_PATH");
53 
54  if( !sp.find_file(fTemplateFile, fROOTfile) )
55  throw cet::exception("Chi2ParticleID") << "cannot find the root template file: \n"
56  << fTemplateFile
57  << "\n bail ungracefully.\n";
58  TFile *file = TFile::Open(fROOTfile.c_str());
59  dedx_range_pro = (TProfile*)file->Get("dedx_range_pro");
60  dedx_range_ka = (TProfile*)file->Get("dedx_range_ka");
61  dedx_range_pi = (TProfile*)file->Get("dedx_range_pi");
62  dedx_range_mu = (TProfile*)file->Get("dedx_range_mu");
63 
64 // std::cout<<"Chi2PIDAlg configuration:"<<std::endl;
65 // std::cout<<"Template file: "<<fROOTfile<<std::endl;
66 // std::cout<<"fUseMedian: "<<fUseMedian<<std::endl;
67 
68  return;
69 }
pid::Chi2PIDAlg::dedx_range_mu
TProfile * dedx_range_mu
muon template
Definition: Chi2PIDAlg.h:50
pid::Chi2PIDAlg::dedx_range_pro
TProfile * dedx_range_pro
proton template
Definition: Chi2PIDAlg.h:47
pid::Chi2PIDAlg::dedx_range_pi
TProfile * dedx_range_pi
pion template
Definition: Chi2PIDAlg.h:49
pid::Chi2PIDAlg::fROOTfile
std::string fROOTfile
Definition: Chi2PIDAlg.h:45
pid::Chi2PIDAlg::fTemplateFile
std::string fTemplateFile
Definition: Chi2PIDAlg.h:42
file
TFile * file
Definition: comparison_ascii.C:5
pid::Chi2PIDAlg::dedx_range_ka
TProfile * dedx_range_ka
kaon template
Definition: Chi2PIDAlg.h:48
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Member Data Documentation

TProfile* pid::Chi2PIDAlg::dedx_range_ka
private

kaon template

Definition at line 48 of file Chi2PIDAlg.h.

Referenced by DoParticleID(), and reconfigure().

TProfile* pid::Chi2PIDAlg::dedx_range_mu
private

muon template

Definition at line 50 of file Chi2PIDAlg.h.

Referenced by DoParticleID(), and reconfigure().

TProfile* pid::Chi2PIDAlg::dedx_range_pi
private

pion template

Definition at line 49 of file Chi2PIDAlg.h.

Referenced by DoParticleID(), and reconfigure().

TProfile* pid::Chi2PIDAlg::dedx_range_pro
private

proton template

Definition at line 47 of file Chi2PIDAlg.h.

Referenced by DoParticleID(), and reconfigure().

std::string pid::Chi2PIDAlg::fROOTfile
private

Definition at line 45 of file Chi2PIDAlg.h.

Referenced by reconfigure().

std::string pid::Chi2PIDAlg::fTemplateFile
private

Definition at line 42 of file Chi2PIDAlg.h.

Referenced by reconfigure().

bool pid::Chi2PIDAlg::fUseMedian
private

Definition at line 43 of file Chi2PIDAlg.h.

Referenced by DoParticleID(), and reconfigure().

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