LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
Chi2PIDAlg.cxx
Go to the documentation of this file.
1 //
3 // Chi2PIDAlg class
4 //
5 // tjyang@fnal.gov
6 //
8 
9 extern "C" {
10 #include <sys/types.h>
11 #include <sys/stat.h>
12 }
13 
14 //#include "RecoBase/Track.h"
19 
20 // ROOT includes
21 #include "TFile.h"
22 #include "TProfile.h"
23 #include "TMath.h"
24 
25 // Framework includes
33 
34 //------------------------------------------------------------------------------
36 {
37  this->reconfigure(pset);
38 }
39 
40 //------------------------------------------------------------------------------
42 {
43 }
44 
45 //------------------------------------------------------------------------------
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 }
70 
71 
72 //------------------------------------------------------------------------------
74  anab::ParticleID &pidOut){
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<double> trkdedx = calo->dEdx();
85  std::vector<double> trkres = calo->ResidualRange();
86  std::vector<double> 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 }
212 
Chi2PIDAlg(fhicl::ParameterSet const &pset)
Definition: Chi2PIDAlg.cxx:35
const std::vector< double > & DeadWireResRC() const
Definition: Calorimetry.h:91
TProfile * dedx_range_mu
muon template
Definition: Chi2PIDAlg.h:50
double fMissingEavg
missing energy from dead wires using average dEdx
Definition: ParticleID.h:34
double fMinChi2
Minimum reduced chi2.
Definition: ParticleID.h:27
const geo::PlaneID & PlaneID() const
Definition: Calorimetry.h:102
TProfile * dedx_range_pro
proton template
Definition: Chi2PIDAlg.h:47
TProfile * dedx_range_pi
pion template
Definition: Chi2PIDAlg.h:49
double fChi2Muon
reduced chi2 using muon template
Definition: ParticleID.h:32
int fNdf
ndf for chi2 test
Definition: ParticleID.h:26
void reconfigure(fhicl::ParameterSet const &pset)
Definition: Chi2PIDAlg.cxx:46
double fDeltaChi2
difference between two lowest reduced chi2&#39;s
Definition: ParticleID.h:28
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::string fROOTfile
Definition: Chi2PIDAlg.h:45
std::string fTemplateFile
Definition: Chi2PIDAlg.h:42
virtual ~Chi2PIDAlg()
Definition: Chi2PIDAlg.cxx:41
float bin[41]
Definition: plottest35.C:14
double TrkPitchC() const
Definition: Calorimetry.h:95
const std::vector< double > & ResidualRange() const
Definition: Calorimetry.h:90
double fChi2Kaon
reduced chi2 using kaon template
Definition: ParticleID.h:30
int fPdg
determined particle ID
Definition: ParticleID.h:25
Utility object to perform functions of association.
double fChi2Proton
reduced chi2 using proton template
Definition: ParticleID.h:29
TFile * file
double fPIDA
PID developed by Bruce Baller.
Definition: ParticleID.h:35
geo::PlaneID fPlaneID
Definition: ParticleID.h:36
TProfile * dedx_range_ka
kaon template
Definition: Chi2PIDAlg.h:48
Float_t e
Definition: plot.C:34
const std::vector< double > & dEdx() const
Definition: Calorimetry.h:88
double fChi2Pion
reduced chi2 using pion template
Definition: ParticleID.h:31
Definition: fwd.h:25
calorimetry
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void DoParticleID(art::Ptr< anab::Calorimetry > calo, anab::ParticleID &pidOut)
Definition: Chi2PIDAlg.cxx:73
double fMissingE
missing energy from dead wires for contained particle
Definition: ParticleID.h:33