13 using namespace bpfit;
33 std::vector<TVector3> traj = track->
Trajectory();
49 for(
unsigned int i = 0;
i+1 < traj.size(); ++
i) {
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);
60 std::vector<geo::CellOnLine> Xhits;
61 std::vector<geo::CellOnLine> Yhits;
71 for(
unsigned int j = 0;
j < Xhits.size(); ++
j) {
75 for(
unsigned int k = 0; k <
chan.size(); ++k) {
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);
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);
107 s.push_back(stot + d1 + d/2.0);
114 for(
unsigned int j = 0;
j < Yhits.size(); ++
j) {
118 for(
unsigned int k = 0; k <
chan.size(); ++k) {
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);
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);
150 s.push_back(stot + d1 + d/2.0);
158 if(
chan.size() !=
dx.size() ||
chan.size() !=
s.size()) abort();
164 for(
unsigned int i = 0;
i < traj.size(); ++
i) {
167 xyz[0] = traj[
i].X();
168 xyz[1] = traj[
i].Y();
169 xyz[2] = traj[
i].Z();
170 path.SetPoint(
i,xyz);
184 for(
unsigned int i = 0;
i <
chan.size(); ++
i) {
194 unsigned int pos = 0;
197 for(
unsigned int j = 0;
j < hits.
size(); ++
j) {
200 if(
chan[
i].Cell() == hits[
j]->Cell() &&
chan[
i].Plane() == hits[
j]->Plane()) {
213 Ctemp = track->
Cell(pos);
214 Wtemp = track->
W(chit);
218 if(rhit.IsCalibrated()) {
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);
225 std::cerr <<
"\n\n\nERROR: PDG code " << pdg <<
" is not supported." 234 dE.push_back(dEtemp);
236 bg.push_back(Btemp*Gtemp);
246 std::vector<double> &Wtemp, std::vector<double> &bgtemp,
247 std::vector<double> &stemp,
double dxTol)
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]);
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]);
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]);
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]);
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]);
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]);
341 for(
unsigned int i = 0;
i <
chan.size(); ++
i) {
342 if(
dE[
i] > 0.0 &&
dx[
i] > dxTol) {
347 return returnedCells;
std::vector< double > bg
value of beta*gamma for the tracked particle at s
Float_t y1[n_points_granero]
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Float_t x1[n_points_granero]
void getBG(std::vector< double > &bgtemp, double dxTol)
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
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.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Calculates dE/dx in cells passed through by a track.
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 push_back(Ptr< U > const &p)
Track parameters (s, X0, KE, beta, ...) along a particular path in the detector.
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
art::ServiceHandle< calib::Calibrator > fCal
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
dEdxCalculator()
Constructor...
A rawdata::RawDigit with channel information decoded.
std::vector< double > W
value of track W for the cell hit
art::ServiceHandle< geo::Geometry > fGeom
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
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 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)