Public Member Functions | Private Attributes | List of all members
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 28 of file Chi2PIDAlg.h.

Constructor & Destructor Documentation

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

Definition at line 26 of file Chi2PIDAlg.cxx.

27 {
28  fTemplateFile = pset.get< std::string >("TemplateFile");
29  fUseMedian = pset.get< bool >("UseMedian");
30  //fCalorimetryModuleLabel = pset.get< std::string >("CalorimetryModuleLabel");
31 
32  cet::search_path sp("FW_SEARCH_PATH");
33 
34  if( !sp.find_file(fTemplateFile, fROOTfile) )
35  throw cet::exception("Chi2ParticleID") << "cannot find the root template file: \n"
36  << fTemplateFile
37  << "\n bail ungracefully.\n";
38  TFile *file = TFile::Open(fROOTfile.c_str());
39  dedx_range_pro = (TProfile*)file->Get("dedx_range_pro");
40  dedx_range_ka = (TProfile*)file->Get("dedx_range_ka");
41  dedx_range_pi = (TProfile*)file->Get("dedx_range_pi");
42  dedx_range_mu = (TProfile*)file->Get("dedx_range_mu");
43 
44 // std::cout<<"Chi2PIDAlg configuration:"<<std::endl;
45 // std::cout<<"Template file: "<<fROOTfile<<std::endl;
46 // std::cout<<"fUseMedian: "<<fUseMedian<<std::endl;
47 }
pid::Chi2PIDAlg::dedx_range_mu
TProfile * dedx_range_mu
muon template
Definition: Chi2PIDAlg.h:51
string
std::string string
Definition: nybbler.cc:12
pid::Chi2PIDAlg::dedx_range_pro
TProfile * dedx_range_pro
proton template
Definition: Chi2PIDAlg.h:48
pid::Chi2PIDAlg::dedx_range_pi
TProfile * dedx_range_pi
pion template
Definition: Chi2PIDAlg.h:50
cache_state.file
file
Definition: cache_state.py:421
pid::Chi2PIDAlg::fROOTfile
std::string fROOTfile
Definition: Chi2PIDAlg.h:46
pid::Chi2PIDAlg::fTemplateFile
std::string fTemplateFile
Definition: Chi2PIDAlg.h:43
pid::Chi2PIDAlg::dedx_range_ka
TProfile * dedx_range_ka
kaon template
Definition: Chi2PIDAlg.h:49
cet::search_path
Definition: search_path.h:36
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Member Function Documentation

anab::ParticleID pid::Chi2PIDAlg::DoParticleID ( const std::vector< art::Ptr< anab::Calorimetry >> &  calo)

Definition at line 61 of file Chi2PIDAlg.cxx.

61  {
62 
63  std::vector<anab::sParticleIDAlgScores> AlgScoresVec;
64  geo::PlaneID plid;
65 
66  for (size_t i_calo = 0; i_calo < calos.size(); i_calo++){
67 
68  art::Ptr<anab::Calorimetry> calo = calos.at(i_calo);
69  if(i_calo == 0)
70  plid = calo->PlaneID();
71  else if(plid != calo->PlaneID())
72  throw cet::exception("Chi2PIDAlg") << "PlaneID mismatch: " << plid << ", " << calo->PlaneID();
73  int npt = 0;
74  double chi2pro = 0;
75  double chi2ka = 0;
76  double chi2pi = 0;
77  double chi2mu = 0;
78  double avgdedx = 0;
79  double PIDA = 0; //by Bruce Baller
80  std::vector<double> vpida;
81  std::vector<float> trkdedx = calo->dEdx();
82  std::vector<float> trkres = calo->ResidualRange();
83  std::vector<float> deadwireresrc = calo->DeadWireResRC();
84 
85  int used_trkres = 0;
86  for (unsigned i = 0; i<trkdedx.size(); ++i){//hits
87  //ignore the first and the last point
88  if (i==0 || i==trkdedx.size()-1) continue;
89  avgdedx += trkdedx[i];
90  if(trkres[i] < 30) {
91  PIDA += trkdedx[i]*pow(trkres[i],0.42);
92  vpida.push_back(trkdedx[i]*pow(trkres[i],0.42));
93  used_trkres++;
94  }
95  if (trkdedx[i]>1000) continue; //protect against large pulse height
96  int bin = dedx_range_pro->FindBin(trkres[i]);
97  if (bin>=1&&bin<=dedx_range_pro->GetNbinsX()){
98  double bincpro = dedx_range_pro->GetBinContent(bin);
99  if (bincpro<1e-6){//for 0 bin content, using neighboring bins
100  bincpro = (dedx_range_pro->GetBinContent(bin-1)+dedx_range_pro->GetBinContent(bin+1))/2;
101  }
102  double bincka = dedx_range_ka->GetBinContent(bin);
103  if (bincka<1e-6){
104  bincka = (dedx_range_ka->GetBinContent(bin-1)+dedx_range_ka->GetBinContent(bin+1))/2;
105  }
106  double bincpi = dedx_range_pi->GetBinContent(bin);
107  if (bincpi<1e-6){
108  bincpi = (dedx_range_pi->GetBinContent(bin-1)+dedx_range_pi->GetBinContent(bin+1))/2;
109  }
110  double bincmu = dedx_range_mu->GetBinContent(bin);
111  if (bincmu<1e-6){
112  bincmu = (dedx_range_mu->GetBinContent(bin-1)+dedx_range_mu->GetBinContent(bin+1))/2;
113  }
114  double binepro = dedx_range_pro->GetBinError(bin);
115  if (binepro<1e-6){
116  binepro = (dedx_range_pro->GetBinError(bin-1)+dedx_range_pro->GetBinError(bin+1))/2;
117  }
118  double bineka = dedx_range_ka->GetBinError(bin);
119  if (bineka<1e-6){
120  bineka = (dedx_range_ka->GetBinError(bin-1)+dedx_range_ka->GetBinError(bin+1))/2;
121  }
122  double binepi = dedx_range_pi->GetBinError(bin);
123  if (binepi<1e-6){
124  binepi = (dedx_range_pi->GetBinError(bin-1)+dedx_range_pi->GetBinError(bin+1))/2;
125  }
126  double binemu = dedx_range_mu->GetBinError(bin);
127  if (binemu<1e-6){
128  binemu = (dedx_range_mu->GetBinError(bin-1)+dedx_range_mu->GetBinError(bin+1))/2;
129  }
130  //double errke = 0.05*trkdedx[i]; //5% KE resolution
131  double errdedx = 0.04231+0.0001783*trkdedx[i]*trkdedx[i]; //resolution on dE/dx
132  errdedx *= trkdedx[i];
133  chi2pro += pow((trkdedx[i]-bincpro)/std::sqrt(pow(binepro,2)+pow(errdedx,2)),2);
134  chi2ka += pow((trkdedx[i]-bincka)/std::sqrt(pow(bineka,2)+pow(errdedx,2)),2);
135  chi2pi += pow((trkdedx[i]-bincpi)/std::sqrt(pow(binepi,2)+pow(errdedx,2)),2);
136  chi2mu += pow((trkdedx[i]-bincmu)/std::sqrt(pow(binemu,2)+pow(errdedx,2)),2);
137  //std::cout<<i<<" "<<trkdedx[i]<<" "<<trkres[i]<<" "<<bincpro<<std::endl;
138  ++npt;
139  }
140  }
141 
142  anab::sParticleIDAlgScores chi2proton;
143  anab::sParticleIDAlgScores chi2kaon;
144  anab::sParticleIDAlgScores chi2pion;
145  anab::sParticleIDAlgScores chi2muon;
146  anab::sParticleIDAlgScores pida_mean;
147  anab::sParticleIDAlgScores pida_median;
148 
149  //anab::ParticleID pidOut;
150  if (npt){
151 
152  chi2proton.fAlgName = "Chi2";
153  chi2proton.fVariableType = anab::kGOF;
154  chi2proton.fTrackDir = anab::kForward;
155  chi2proton.fAssumedPdg = 2212;
156  chi2proton.fPlaneMask = GetBitset(calo->PlaneID());
157  chi2proton.fNdf = npt;
158  chi2proton.fValue = chi2pro/npt;
159 
160  chi2muon.fAlgName = "Chi2";
161  chi2muon.fVariableType = anab::kGOF;
162  chi2muon.fTrackDir = anab::kForward;
163  chi2muon.fAssumedPdg = 13;
164  chi2muon.fPlaneMask = GetBitset(calo->PlaneID());
165  chi2muon.fNdf = npt;
166  chi2muon.fValue = chi2mu/npt;
167 
168  chi2kaon.fAlgName = "Chi2";
169  chi2kaon.fVariableType = anab::kGOF;
170  chi2kaon.fTrackDir = anab::kForward;
171  chi2kaon.fAssumedPdg = 321;
172  chi2kaon.fPlaneMask = GetBitset(calo->PlaneID());
173  chi2kaon.fNdf = npt;
174  chi2kaon.fValue = chi2ka/npt;
175 
176  chi2pion.fAlgName = "Chi2";
177  chi2pion.fVariableType = anab::kGOF;
178  chi2pion.fTrackDir = anab::kForward;
179  chi2pion.fAssumedPdg = 211;
180  chi2pion.fPlaneMask = GetBitset(calo->PlaneID());
181  chi2pion.fNdf = npt;
182  chi2pion.fValue = chi2pi/npt;
183 
184  AlgScoresVec.push_back(chi2proton);
185  AlgScoresVec.push_back(chi2muon);
186  AlgScoresVec.push_back(chi2kaon);
187  AlgScoresVec.push_back(chi2pion);
188  }
189 
190  //if (trkdedx.size()) pidOut.fPIDA = PIDA/trkdedx.size();
191  if(used_trkres > 0){
192  if (fUseMedian){
193  pida_median.fAlgName = "PIDA_median";
194  pida_median.fVariableType = anab::kPIDA;
195  pida_median.fTrackDir = anab::kForward;
196  pida_median.fValue = TMath::Median(vpida.size(), &vpida[0]);
197  pida_median.fPlaneMask = GetBitset(calo->PlaneID());
198  AlgScoresVec.push_back(pida_median);
199  }
200  else{ // use mean
201  pida_mean.fAlgName = "PIDA_mean";
202  pida_mean.fVariableType = anab::kPIDA;
203  pida_mean.fTrackDir = anab::kForward;
204  pida_mean.fValue = PIDA/used_trkres;
205  pida_mean.fPlaneMask = GetBitset(calo->PlaneID());
206  AlgScoresVec.push_back(pida_mean);
207  }
208  }
209 
210  }
211 
212  anab::ParticleID pidOut(AlgScoresVec, plid);
213 
214  return pidOut;
215 
216 }
pid::Chi2PIDAlg::dedx_range_mu
TProfile * dedx_range_mu
muon template
Definition: Chi2PIDAlg.h:51
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
anab::Calorimetry::PlaneID
const geo::PlaneID & PlaneID() const
Definition: Calorimetry.h:116
pid::Chi2PIDAlg::dedx_range_pro
TProfile * dedx_range_pro
proton template
Definition: Chi2PIDAlg.h:48
pid::Chi2PIDAlg::dedx_range_pi
TProfile * dedx_range_pi
pion template
Definition: Chi2PIDAlg.h:50
anab::sParticleIDAlgScores::fValue
float fValue
Result of Particle ID algorithm/test.
Definition: ParticleID.h:28
anab::Calorimetry::DeadWireResRC
const std::vector< float > & DeadWireResRC() const
Definition: Calorimetry.h:104
anab::sParticleIDAlgScores::fPlaneMask
std::bitset< 8 > fPlaneMask
Bitset for PlaneID used by algorithm, allowing for multiple planes and up to 8 total planes...
Definition: ParticleID.h:29
anab::Calorimetry::ResidualRange
const std::vector< float > & ResidualRange() const
Definition: Calorimetry.h:103
anab::sParticleIDAlgScores::fAlgName
std::string fAlgName
< determined particle ID
Definition: ParticleID.h:23
anab::sParticleIDAlgScores::fNdf
int fNdf
Number of degrees of freedom used by algorithm, if applicable. Set to -9999 by default.
Definition: ParticleID.h:26
anab::sParticleIDAlgScores::fTrackDir
kTrackDir fTrackDir
Track direction enum: defined in ParticleID_VariableTypeEnums.h. Set to kNoDirection by default...
Definition: ParticleID.h:25
anab::sParticleIDAlgScores::fVariableType
kVariableType fVariableType
Variable type enum: defined in ParticleID_VariableTypeEnums.h. Set to kNotSet by default.
Definition: ParticleID.h:24
e
const double e
Definition: gUpMuFluxGen.cxx:165
pid::Chi2PIDAlg::GetBitset
std::bitset< 8 > GetBitset(geo::PlaneID planeID)
Definition: Chi2PIDAlg.cxx:50
anab::Calorimetry::dEdx
const std::vector< float > & dEdx() const
Definition: Calorimetry.h:101
bin
QTextStream & bin(QTextStream &s)
Definition: qtextstream.cpp:2006
pid::Chi2PIDAlg::dedx_range_ka
TProfile * dedx_range_ka
kaon template
Definition: Chi2PIDAlg.h:49
art::Ptr
Definition: fwd.h:31
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
anab::sParticleIDAlgScores::fAssumedPdg
int fAssumedPdg
PDG of particle hypothesis assumed by algorithm, if applicable. Set to 0 by default.
Definition: ParticleID.h:27
std::bitset< 8 > pid::Chi2PIDAlg::GetBitset ( geo::PlaneID  planeID)

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

Definition at line 50 of file Chi2PIDAlg.cxx.

50  {
51 
52  std::bitset<8> thisBitset;
53 
54  thisBitset.set(planeID.Plane);
55 
56  return thisBitset;
57 
58 }
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493

Member Data Documentation

TProfile* pid::Chi2PIDAlg::dedx_range_ka
private

kaon template

Definition at line 49 of file Chi2PIDAlg.h.

TProfile* pid::Chi2PIDAlg::dedx_range_mu
private

muon template

Definition at line 51 of file Chi2PIDAlg.h.

TProfile* pid::Chi2PIDAlg::dedx_range_pi
private

pion template

Definition at line 50 of file Chi2PIDAlg.h.

TProfile* pid::Chi2PIDAlg::dedx_range_pro
private

proton template

Definition at line 48 of file Chi2PIDAlg.h.

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

Definition at line 46 of file Chi2PIDAlg.h.

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

Definition at line 43 of file Chi2PIDAlg.h.

bool pid::Chi2PIDAlg::fUseMedian
private

Definition at line 44 of file Chi2PIDAlg.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11