dEdxCalculator.cxx
Go to the documentation of this file.
1 ///
2 /// \file dEdxCalculator.cxx
3 /// \brief Calculates dE/dx in cells passed through by a track.
4 ///
5 
6 #include "RecoBase/CellHit.h"
7 #include "RecoBase/RecoHit.h"
8 #include "RecoBase/Track.h"
9 
12 
13 using namespace bpfit;
14 
15 //static const double gsMmuon = 105.658E-3;
16 //static const double gsMpion = 139.570E-3;
17 //static const double gsMproton = 938.272E-3;
18 
19 //......................................................................
20 
22 
23 //......................................................................
24 
26 {
27 
28  // Reset.
29  this->clearVectors();
30 
31  // get the cells in the track and get the trajectory points
33  std::vector<TVector3> traj = track->Trajectory();
34 
35  // keep track of trajectory length so far
36  double stot = 0.0;
37 
38  //
39  // For each track, we will follow along the trajectory segments and make
40  // a list of all of the cells that the track passed through. For each of
41  // these cells, we will compute dx as the length between the entry point
42  // and the exit point. Since it is possible for a trajectory segment to
43  // end inside a cell (with the next segment then continuing on to exit
44  // the cell) we will sum the computed dx's for cells that appear on our
45  // list more than once.
46  //
47 
48  // Loop over trajectory points:
49  for(unsigned int i = 0; i+1 < traj.size(); ++i) {
50 
51  // add length of previous segment to total path length
52  if(i > 0) {
53  double x0 = traj[i].X() - traj[i-1].X();
54  double y0 = traj[i].Y() - traj[i-1].Y();
55  double z0 = traj[i].Z() - traj[i-1].Z();
56  double d0 = sqrt(x0*x0 + y0*y0 + z0*z0);
57  stot += d0;
58  }
59 
60  std::vector<geo::CellOnLine> Xhits;
61  std::vector<geo::CellOnLine> Yhits;
62 
63  // make list of cells passed through by this segment of the trajectory
64  fGeom->CountCellsOnLine(traj[i],traj[i+1],Xhits,Yhits);
65 
66  // Check if hits are already in my list. If it is not on the list,
67  // add it and its info. If it is, then only add the additional
68  // computed dx.
69  bool onlist = false;
70  unsigned int pos = 0;
71  for(unsigned int j = 0; j < Xhits.size(); ++j) {
72 
73  onlist = false;
74  pos = 0;
75  for(unsigned int k = 0; k < chan.size(); ++k) {
76  if(Xhits[j].chan == chan[k]) {
77  onlist = true;
78  pos = k;
79  break;
80  }
81  } // end loop over elements of chan
82 
83  // compute path length through the cell (for dx)
84  double x = Xhits[j].entry.X() - Xhits[j].exit.X();
85  double y = Xhits[j].entry.Y() - Xhits[j].exit.Y();
86  double z = Xhits[j].entry.Z() - Xhits[j].exit.Z();
87  double d = sqrt(x*x + y*y + z*z);
88 
89  // compute length from trajectory point i to entry point of cell
90  // (for proper calculation of path length s for this cell)
91  double x1 = Xhits[j].entry.X() - traj[i].X();
92  double y1 = Xhits[j].entry.Y() - traj[i].Y();
93  double z1 = Xhits[j].entry.Z() - traj[i].Z();
94  double d1 = sqrt(x1*x1 + y1*y1 + z1*z1);
95 
96  if(onlist) {
97  dx [pos] += d;
98  s [pos] += d/2.0;
99  }
100  else {
101  chan.push_back(Xhits[j].chan);
102  dx.push_back(d);
103  // The path length for this hit is taken to be the trajectory length
104  // (up to the i-th trajectory point) plus the distance from the i-th
105  // trajectory point to the entrance to the cell, plus half the distance
106  // travelled through the cell.
107  s.push_back(stot + d1 + d/2.0);
108  }
109 
110  } // end loop over elements of Xhits
111 
112  onlist = false;
113  pos = 0;
114  for(unsigned int j = 0; j < Yhits.size(); ++j) {
115 
116  onlist = false;
117  pos = 0;
118  for(unsigned int k = 0; k < chan.size(); ++k) {
119  if(Yhits[j].chan == chan[k]) {
120  onlist = true;
121  pos = k;
122  break;
123  }
124  } // end loop over elements of chan
125 
126  // compute path length through the cell
127  double x = Yhits[j].entry.X() - Yhits[j].exit.X();
128  double y = Yhits[j].entry.Y() - Yhits[j].exit.Y();
129  double z = Yhits[j].entry.Z() - Yhits[j].exit.Z();
130  double d = sqrt(x*x + y*y + z*z);
131 
132  // compute length from trajectory point i to entry point of cell
133  // (for proper calculation of path length s for this cell)
134  double x1 = Yhits[j].entry.X() - traj[i].X();
135  double y1 = Yhits[j].entry.Y() - traj[i].Y();
136  double z1 = Yhits[j].entry.Z() - traj[i].Z();
137  double d1 = sqrt(x1*x1 + y1*y1 + z1*z1);
138 
139  if(onlist) {
140  dx[pos] += d;
141  s [pos] += d/2.0;
142  }
143  else {
144  chan.push_back(Yhits[j].chan);
145  dx.push_back(d);
146  // The path length for this hit is taken to be the trajectory length
147  // (up to the i-th trajectory point) plus the distance from the i-th
148  // trajectory point to the entrance to the cell, plus half the distance
149  // travelled through the cell.
150  s.push_back(stot + d1 + d/2.0);
151  }
152 
153  } // end loop over elements of Yhits
154 
155  } // end loop over trajectory points (i)
156 
157  // check that we haven't messed up
158  if(chan.size() != dx.size() || chan.size() != s.size()) abort();
159 
160 
161 
162  // make the bpfit::path object to compute momentum along the trajectory
163  bpfit::Path path(*fGeom, traj.size());
164  for(unsigned int i = 0; i < traj.size(); ++i) {
165  // set the points in path using the trajectory points
166  double xyz[3];
167  xyz[0] = traj[i].X();
168  xyz[1] = traj[i].Y();
169  xyz[2] = traj[i].Z();
170  path.SetPoint(i,xyz);
171  }
172 
173 
174 
175 
176  //
177  // Now populate the list of energies and momenta. Since a trajectory can pass
178  // through a cell without there being a cell hit, we will loop over all cells
179  // that the trajectory passed through and only compute dE/dx for the ones that
180  // have a corresponding cell hit on the track.
181  //
182 
183  // loop over cells in chan and check for a cell hit in the list of hits in the track
184  for(unsigned int i = 0; i < chan.size(); ++i) {
185 
186  art::Ptr<rb::CellHit> Ctemp = track->Cell(0);
187  double dEtemp = -1;
188  double Wtemp = -1e3;
189  double ptemp = -1;
190  double Btemp = -1;
191  double Gtemp = 1;
192  double Ttemp = -1;
193  bool match = false;
194  unsigned int pos = 0;
195 
196  // look for this cell in the list of hits in the track
197  for(unsigned int j = 0; j < hits.size(); ++j) {
198 
199  // match by plane and cell number (can this assumption lead to trouble???)
200  if(chan[i].Cell() == hits[j]->Cell() && chan[i].Plane() == hits[j]->Plane()) {
201  match = true;
202  pos = j;
203  break;
204  }
205 
206  } // end loop over hits on track (j)
207 
208  // A match was found, compute the real energy of the cell correcting for attenuation calibration.
209  if(match) {
210 
211  // make reco hit from the cell hit
212  const rb::CellHit *chit = track->Cell(pos).get();
213  Ctemp = track->Cell(pos);
214  Wtemp = track->W(chit);
215  rb::RecoHit rhit(fCal->MakeRecoHit(*chit,Wtemp));
216 
217  // if reco hit is calibrated, get GeV to use as dE and compute the momentum using path
218  if(rhit.IsCalibrated()) {
219  dEtemp = rhit.GeV();
220 
221  if (pdg == 13) path.MuonParams(s[i], &Ttemp, &ptemp, &Btemp, &Gtemp);
222  else if(pdg == 211) path.PionParams(s[i], &Ttemp, &ptemp, &Btemp, &Gtemp);
223  else if(pdg == 2212) path.ProtonParams(s[i], &Ttemp, &ptemp, &Btemp, &Gtemp);
224  else {
225  std::cerr << "\n\n\nERROR: PDG code " << pdg << " is not supported."
226  << "\tABORTING...";
227  abort();
228  }
229 
230  }
231 
232  }
233 
234  dE.push_back(dEtemp);
235  W .push_back(Wtemp);
236  bg.push_back(Btemp*Gtemp);
237  cells.push_back(Ctemp);
238 
239  } // end loop over cells in chan list (i)
240 
241 }
242 
243 //......................................................................
244 
245 void dEdxCalculator::getResults(std::vector<double> &dEtemp, std::vector<double> &dxtemp,
246  std::vector<double> &Wtemp, std::vector<double> &bgtemp,
247  std::vector<double> &stemp, double dxTol)
248 {
249 
250  // Construct the list of successfully calculated dE/dx values.
251  for(unsigned int i = 0; i < chan.size(); ++i) {
252  if(dE[i] > 0.0 && dx[i] > dxTol) {
253  dEtemp.push_back(dE[i]);
254  dxtemp.push_back(dx[i]);
255  Wtemp .push_back(W[i]);
256  bgtemp.push_back(bg[i]);
257  stemp .push_back(s[i]);
258  }
259  }
260 
261 }
262 
263 //......................................................................
264 
265 void dEdxCalculator::getDE(std::vector<double> &dEtemp, double dxTol)
266 {
267 
268  // Construct the list of successfully calculated values.
269  for(unsigned int i = 0; i < chan.size(); ++i) {
270  if(dE[i] > 0.0 && dx[i] > dxTol) {
271  dEtemp.push_back(dE[i]);
272  }
273  }
274 
275 }
276 
277 //......................................................................
278 
279 void dEdxCalculator::getDX(std::vector<double> &dxtemp, double dxTol)
280 {
281 
282  // Construct the list of successfully calculated values.
283  for(unsigned int i = 0; i < chan.size(); ++i) {
284  if(dE[i] > 0.0 && dx[i] > dxTol) {
285  dxtemp.push_back(dx[i]);
286  }
287  }
288 
289 }
290 
291 //......................................................................
292 
293 void dEdxCalculator::getW(std::vector<double> &Wtemp, double dxTol)
294 {
295 
296  // Construct the list of successfully calculated values.
297  for(unsigned int i = 0; i < chan.size(); ++i) {
298  if(dE[i] > 0.0 && dx[i] > dxTol) {
299  Wtemp.push_back(W[i]);
300  }
301  }
302 
303 }
304 
305 //......................................................................
306 
307 void dEdxCalculator::getBG(std::vector<double> &bgtemp, double dxTol)
308 {
309 
310  // Construct the list of successfully calculated values.
311  for(unsigned int i = 0; i < chan.size(); ++i) {
312  if(dE[i] > 0.0 && dx[i] > dxTol) {
313  bgtemp.push_back(bg[i]);
314  }
315  }
316 
317 }
318 
319 //......................................................................
320 
321 void dEdxCalculator::getS(std::vector<double> &stemp, double dxTol)
322 {
323 
324  // Construct the list of successfully calculated values.
325  for(unsigned int i = 0; i < chan.size(); ++i) {
326  if(dE[i] > 0.0 && dx[i] > dxTol) {
327  stemp.push_back(s[i]);
328  }
329  }
330 
331 }
332 
333 //......................................................................
334 
336 {
337 
338  art::PtrVector<rb::CellHit> returnedCells;
339 
340  // Construct the list of cells used to calculate the dE/dx values.
341  for(unsigned int i = 0; i < chan.size(); ++i) {
342  if(dE[i] > 0.0 && dx[i] > dxTol) {
343  returnedCells.push_back(cells[i]);
344  }
345  }
346 
347  return returnedCells;
348 
349 }
350 
351 //......................................................................
352 
354 {
355 
356  // Clear all vectors used.
357  chan .clear();
358  cells.clear();
359  dE .clear();
360  dx .clear();
361  W .clear();
362  s .clear();
363  bg .clear();
364 
365 }
366 
367 ////////////////////////////////////////////////////////////////////////
std::vector< double > bg
value of beta*gamma for the tracked particle at s
Float_t y1[n_points_granero]
Definition: compare.C:5
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Track.cxx:281
Float_t x1[n_points_granero]
Definition: compare.C:5
void getBG(std::vector< double > &bgtemp, double dxTol)
T sqrt(T number)
Definition: d0nt_math.hpp:156
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
Definition: event.h:19
OStream cerr
Definition: OStream.cxx:7
void getS(std::vector< double > &stemp, double dxTol)
void computeDEDX(art::Ptr< rb::Track > &track, int pdg=13)
Compute dE/dx for all cells along this track and the fitsum object that went with it...
std::vector< double > dE
computed value of energy for the cell hit
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
Calculates dE/dx in cells passed through by a track.
const XML_Char * s
Definition: expat.h:262
void clearVectors()
Clear all internally used vectors.
void CountCellsOnLine(double X1, double Y1, double Z1, double X2, double Y2, double Z2, std::vector< OfflineChan > &Xhitsonline, std::vector< OfflineChan > &Yhitsonline)
void hits()
Definition: readHits.C:15
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
Track parameters (s, X0, KE, beta, ...) along a particular path in the detector.
Definition: Path.h:18
std::vector< double > dx
computed value of the track path length through the cell hit
art::PtrVector< rb::CellHit > cells
list of cell hits on the track
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
art::ServiceHandle< calib::Calibrator > fCal
z
Definition: test.py:28
size_type size() const
Definition: PtrVector.h:308
const std::string path
Definition: plot_BEN.C:43
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
dEdxCalculator()
Constructor...
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
float GeV() const
Definition: RecoHit.cxx:69
std::vector< double > W
value of track W for the cell hit
art::ServiceHandle< geo::Geometry > fGeom
T const * get() const
Definition: Ptr.h:321
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
Definition: Track.h:86
std::vector< geo::OfflineChan > chan
offline channel associated with a cell hit on a track
void getDE(std::vector< double > &dEtemp, double dxTol)
...can also get the individual results...
void getW(std::vector< double > &Wtemp, double dxTol)
art::PtrVector< rb::CellHit > getCellHits(double dxTol)
Return the list of rb::CellHits used in the dE/dx calculation.
void clear()
Definition: PtrVector.h:537
void getResults(std::vector< double > &dEtemp, std::vector< double > &dxtemp, std::vector< double > &Wtemp, std::vector< double > &bgtemp, std::vector< double > &stemp, double dxTol)
Return the results of the dE/dx calculation.
void getDX(std::vector< double > &dxtemp, double dxTol)