ReMId.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file ReMId.cxx
3 // \brief A pid object for muon identification
4 // \version $Id: ReMId.cxx,v 1.1.1.1 2012-09-26 18:38:44 raddatz Exp $
5 // \author Nicholas Raddatz - raddatz@physics.umn.edu
6 ////////////////////////////////////////////////////////////////////////
7 #include <limits>
8 #include <iostream>
9 
10 #include "ReMId.h"
11 #include "TVector3.h"
12 
14 #include "cetlib_except/exception.h"
15 
16 namespace remid
17 {
18 
19  //......................................
21  : rb::PID()
22  {
23  SetMuonDedxLL(0);
24  SetPionDedxLL(0);
25  SetMuonScatLL(0);
26  SetMuonScatLL(0);
27  }
28 
29  //......................................
30  ReMId::ReMId(int pdg, double val)
31  : rb::PID(pdg,val)
32  {
33  SetMuonDedxLL(0);
34  SetPionDedxLL(0);
35  SetMuonScatLL(0);
36  SetMuonScatLL(0);
37  }
38 
39  //......................................
40  void ReMId::AddScatter(TVector3 pos,double cosScat)
41  {
42  fscatLocation.push_back(pos);
43  fcosScatAngle.push_back(cosScat);
44  }
45 
46  //.....................................
47  void ReMId::AddDedx(double dedx, double distFromStart, bool used, bool vert)
48  {
49  fdedx.push_back(dedx);
50  ffracDistFromStart.push_back(distFromStart);
51  fdedxused.push_back(used);
52  fdedxvert.push_back(vert);
53  }
54 
55  //.....................................
56  void ReMId::SetMuonDedxLL(double ll)
57  {
58  fMuonDedxLL = ll;
59  }
60 
61  //.....................................
62  void ReMId::SetPionDedxLL(double ll)
63  {
64  fPionDedxLL = ll;
65  }
66 
67  //.....................................
68  double ReMId::AvgDedx() const
69  {
70  if(fdedx.size() == 0){ return 0;}
71  else{
72  double totDedx = 0;
73  for(unsigned int i = 0; i<fdedx.size();++i){
74  totDedx+=fdedx[i];
75  }
76  return totDedx/((double)fdedx.size());
77  }
78  }
79 
80  //.....................................
81  double ReMId::AvgDedxUsed() const
82  {
83  int ndedxused = 0;
84  double totDedx = 0;
85  for(unsigned int i = 0; i < fdedx.size(); ++i){
86  if(fdedxused[i]){
87  ++ndedxused;
88  totDedx+=fdedx[i];
89  }
90  }
91  if(ndedxused == 0){ return 0;}
92  else{return totDedx/((double)ndedxused); }
93  }
94 
95  //.....................................
96  unsigned int ReMId::NScatters() const
97  {
98  return fscatLocation.size();
99  }
100 
101  //.....................................
102  unsigned int ReMId::NDedx() const
103  {
104  return fdedx.size();
105  }
106 
107  //.....................................
108  TVector3 ReMId::ScatLocation(unsigned int i) const
109  {
110  return fscatLocation[i];
111  }
112 
113  //.....................................
114  double ReMId::CosScat(unsigned int i) const
115  {
116  return fcosScatAngle[i];
117  }
118 
119  //.....................................
120  double ReMId::DedxLocation(unsigned int i) const
121  {
122  return ffracDistFromStart[i];
123  }
124 
125  //.....................................
126  double ReMId::Dedx(unsigned int i) const
127  {
128  return fdedx[i];
129  }
130 
131  //.....................................
132  bool ReMId::DedxUsed(unsigned int i) const
133  {
134  return fdedxused[i];
135  }
136 
137  //.....................................
138  bool ReMId::DedxVertex(unsigned int i) const
139  {
140  return fdedxvert[i];
141  }
142 
143  //.....................................
144  double ReMId::MuonDedxLL() const
145  {
146  return fMuonDedxLL;
147  }
148 
149  //.....................................
150  double ReMId::PionDedxLL() const
151  {
152  return fPionDedxLL;
153  }
154 
155  //.....................................
156  void ReMId::SetMuonScatLL(double ll)
157  {
158  fMuonScatLL = ll;
159  }
160 
161  //.....................................
162  void ReMId::SetPionScatLL(double ll)
163  {
164  fPionScatLL = ll;
165  }
166 
167  //.....................................
168  void ReMId::SetDedxUsed(unsigned int idx, bool used)
169  {
170  fdedxused[idx] = used;
171  }
172 
173  //.....................................
174  void ReMId::SetDedxVertex(unsigned int idx, bool vert)
175  {
176  fdedxvert[idx] = vert;
177  }
178 
179  //.....................................
180  void ReMId::SetContained(bool cont)
181  {
182  fCont = cont;
183  }
184 
185  //.....................................
186  double ReMId::MuonScatLL() const
187  {
188  return fMuonScatLL;
189  }
190 
191  //.....................................
192  double ReMId::PionScatLL() const
193  {
194  return fPionScatLL;
195  }
196 
197  //.....................................
198  /// Return the dE/dx separation variable used as an input to the kNN that determines the pid value
199  double ReMId::DedxSeparation() const
200  {
201  return fMuonDedxLL-fPionDedxLL;
202  }
203 
204  //.....................................
205  /// Return the scattering separation variable used as an input to the kNN that determines the pid value
206  double ReMId::ScatSeparation() const
207  {
208  return fMuonScatLL-fPionScatLL;
209  }
210 
211  //.....................................
212  /// Return the measurement fraction variable used as an input to the kNN that determines the pid value
214  {
215  double frac = -1;
216  if(this->NDedx() > 0){
217  frac = ((double)this->NMeasurementPlanes())/((double)this->NDedx());
218  }
219  return frac;
220  }
221 
222  //.....................................
224  {
225  int npl = 0;
226  for(unsigned int i = 0; i < fdedxused.size(); ++i){
227  if(fdedxused[i]){ ++npl; }
228  }
229  return npl;
230  }
231 
232  //.....................................
234  {
235  int npl = 0;
236  for(unsigned int i = 0; i< fdedxvert.size(); ++i){
237  if(fdedxvert[i]){ ++npl; }
238  }
239  return npl;
240  }
241 
242  //.....................................
243  bool ReMId::Contained() const
244  {
245  return fCont;
246  }
247 
248  //.....................................
249  unsigned int HighestPIDTrack(
250  const std::vector< art::Ptr<rb::Track> >& sliceTracks,
251  const std::string& remidModuleLabel, const art::Event& e)
252  {
253 
254  art::InputTag tag(remidModuleLabel);
255  return remid::HighestPIDTrack(sliceTracks, tag, e);
256  }
257 
258  //.....................................
259  unsigned int HighestPIDTrack(std::vector< art::Ptr<rb::Track> > const& sliceTracks,
260  art::InputTag const& remidTag,
261  art::Event const& e)
262  {
263  int bestTrack = -1;
264  double highestPID = -2.0;
265  double bestTrackLength = 0.0;
266  double bestStartZ = std::numeric_limits<double>::max();
267  double bestADC = std::numeric_limits<double>::max();
268 
269  art::FindOneP<remid::ReMId> tracksToReMId(sliceTracks, e, remidTag);
270 
271  for(size_t c = 0; c < sliceTracks.size(); ++c){
272 
273  geo::View_t view = sliceTracks[c]->View();
274 
275  //Only use 3D tracks
276  if (view != geo::kXorY){
277  continue;
278  } //End of loop over 2d tracks
279 
280  if(tracksToReMId.isValid()){
281  art::Ptr<remid::ReMId> trackReMId = tracksToReMId.at(c);
282 
283  // If the PID is lower, we can bail, nothing to update
284  if(trackReMId->Value() < highestPID) continue;
285 
286  if(trackReMId->Value() == highestPID)
287  { // If the new guy is shorter, we can bail
288  if (sliceTracks[c]->TotalLength() < bestTrackLength) continue;
289 
290  // If equal, need to fall to the next level, low start Z
291  if (sliceTracks[c]->TotalLength() == bestTrackLength){
292 
293  // If start Z is higher, we can bail
294  if (sliceTracks[c]->Start().Z() > bestStartZ) continue;
295 
296  // If equal, need to fall to the next level, total ADC
297  // You might think this is stupid, but it's equivalent to flipping
298  // a coin which always lands on the same side. This is good, this
299  // will happen so rarely that you should just pee up a rope.
300  // Use lower ADC
301  if (sliceTracks[c]->Start().Z() == bestStartZ){
302 
303  if(sliceTracks[c]->TotalADC() > bestADC) continue;
304 
305  if(sliceTracks[c]->TotalADC() == bestADC)
306  throw cet::exception("ReMId Matched Track Exception")<<std::endl
307  << "Can't break a tie between two ReMId/tracks. "
308  << "This should have never happend, but it did. "
309  << "The numu group is going to have to reflect on this "
310  << "and find a way to fix it." << std::endl;
311 
312  }
313 
314 
315  }
316  }
317 
318 
319 
320  highestPID = trackReMId->Value();
321  bestTrackLength = sliceTracks[c]->TotalLength();
322  bestStartZ = sliceTracks[c]->Start().Z();
323  bestADC = sliceTracks[c]->TotalADC();
324  bestTrack = c;
325 
326 
327  }//End of valid ReMId object
328 
329  }//End of for loop over tracks
330 
331 
332  if (bestTrack < 0)
333  return 999;
334 
335  return (unsigned int)bestTrack;
336 
337 
338  }
339 
340 
341 }
double PionDedxLL() const
Definition: ReMId.cxx:150
std::vector< double > fdedx
Definition: ReMId.h:69
double Dedx(unsigned int i) const
Definition: ReMId.cxx:126
void SetPionScatLL(double ll)
Definition: ReMId.cxx:162
unsigned int NScatters() const
Definition: ReMId.cxx:96
double MuonDedxLL() const
Definition: ReMId.cxx:144
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
PID()
Definition: PID.cxx:13
std::vector< int > fdedxvert
Definition: ReMId.h:70
double fMuonDedxLL
Definition: ReMId.h:63
double ScatSeparation() const
Return the scattering separation variable used as an input to the kNN that determines the pid value...
Definition: ReMId.cxx:206
bool fCont
Definition: ReMId.h:65
void SetDedxVertex(unsigned int idx, bool vert)
Definition: ReMId.cxx:174
double MuonScatLL() const
Definition: ReMId.cxx:186
std::vector< int > fdedxused
Definition: ReMId.h:71
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void SetDedxUsed(unsigned int idx, bool used)
Definition: ReMId.cxx:168
std::vector< TVector3 > fscatLocation
Definition: ReMId.h:67
void SetContained(bool cont)
Definition: ReMId.cxx:180
double AvgDedx() const
Definition: ReMId.cxx:68
double Value() const
Definition: PID.h:22
double fPionDedxLL
Definition: ReMId.h:64
bool DedxVertex(unsigned int i) const
Definition: ReMId.cxx:138
void SetMuonDedxLL(double ll)
Definition: ReMId.cxx:56
void AddDedx(double dedx, double distFromStart, bool used=true, bool vert=false)
Definition: ReMId.cxx:47
int NMeasurementPlanes() const
Definition: ReMId.cxx:223
TVector3 ScatLocation(unsigned int i) const
Definition: ReMId.cxx:108
double CosScat(unsigned int i) const
Definition: ReMId.cxx:114
std::vector< double > ffracDistFromStart
Definition: ReMId.h:72
void SetPionDedxLL(double ll)
Definition: ReMId.cxx:62
int NVertexPlanes() const
Definition: ReMId.cxx:233
bool DedxUsed(unsigned int i) const
Definition: ReMId.cxx:132
bool Contained() const
Definition: ReMId.cxx:243
double fMuonScatLL
Definition: ReMId.h:61
double frac(double x)
Fractional part.
void AddScatter(TVector3 pos, double cosScat)
Definition: ReMId.cxx:40
Perform a "2 point" Hough transform on a collection of hits.
void SetMuonScatLL(double ll)
Definition: ReMId.cxx:156
double AvgDedxUsed() const
Definition: ReMId.cxx:81
unsigned int NDedx() const
Definition: ReMId.cxx:102
double fPionScatLL
Definition: ReMId.h:62
double DedxLocation(unsigned int i) const
Definition: ReMId.cxx:120
double MeasurementFraction() const
Return the measurement fraction variable used as an input to the kNN that determines the pid value...
Definition: ReMId.cxx:213
std::vector< double > fcosScatAngle
Definition: ReMId.h:68
double PionScatLL() const
Definition: ReMId.cxx:192
A PID for muons.
Definition: FillPIDs.h:10
Float_t e
Definition: plot.C:35
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
double DedxSeparation() const
Return the dE/dx separation variable used as an input to the kNN that determines the pid value...
Definition: ReMId.cxx:199
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
Definition: ReMId.cxx:249
enum BeamMode string