NDRecoFxs.cxx
Go to the documentation of this file.
1 #include "NDReco/NDRecoFxs.h"
2 
4 #include "Geometry/Geometry.h"
6 #include "GeometryObjects/Geo.h"
9 #include "RecoBase/RecoHit.h"
10 
12 
13 #include "TMath.h"
14 #include "TVector3.h"
15 
16 
17 namespace ndreco
18 {
19 
20  //-------------------------------------------------------
21  NDRecoFxs::NDRecoFxs() //constructor
22  {
23  }
24 
25  //-------------------------------------------------------
26  NDRecoFxs::~NDRecoFxs() //destructor
27  {
28  }
29 
30  //-------------------------------------------------------
31 
32  float NDRecoFxs::getTrackScatt(const art::Ptr<rb::Track> track, TVector3 vert, float dist) {
33  double tscattmax=0, tscatttot=0, tscatt;
34  TVector3 traj1, traj2, dir, olddir;
35 
36  for(unsigned int p=0;p<track->NTrajectoryPoints()-1;p++) {
37  traj1 = track->TrajectoryPoint(p);
38  traj2 = track->TrajectoryPoint(p+1);
39  dir = (traj2-traj1).Unit();
40  if((dir-olddir).Mag()<1e-7) continue;
41  if(p>0&&dir.Mag()<=1.0&&olddir.Mag()<=1.0&&(traj1-vert).Mag()>dist&&(traj2-vert).Mag()>dist) {
42  tscatt = std::acos(dir.Dot(olddir)) * 180 / TMath::Pi();
43  if(tscatt > tscattmax) tscattmax = tscatt;
44  tscatttot += tscatt;
45  }
46  olddir = dir;
47  }
48  return tscatttot; // calculated both maximum and total scatter, total seems more useful, ignore max for now
49  }
50 
51  //-------------------------------------------------------
52  // crude dE/dx: dE is GeV of hits, dx is a track length. Complicated by excluding region within dist cm of vertex
53  float NDRecoFxs::getTrackDedx(const art::Ptr<rb::Track> track, TVector3 vert, float dist) {
55  float gev = 0, length = 0;
56  TVector3 traj1, traj2;
57 
58  art::PtrVector< rb::CellHit > const& trackHits = track->AllCells();
59 
60  rb::Cluster survivingCluster; // make cluster of hits further than dist away from vertex
61 
62  for (unsigned int trackHit = 0; trackHit < trackHits.size(); ++trackHit){ //fill cluster
63  const geo::PlaneGeo* p = geom->Plane(trackHits[trackHit]->Plane());
64  double xyz[3] = {0.,0.,0.}, junk = 0;
65  p->Cell(trackHits[trackHit]->Cell())->GetCenter(xyz, junk);
66  float dist2d = 0;
67  if(trackHits[trackHit]->View()==geo::kX) dist2d = sqrt( pow(xyz[0]-vert.X(),2) + pow(xyz[2]-vert.Z(),2));
68  else if(trackHits[trackHit]->View()==geo::kY) dist2d = sqrt( pow(xyz[1]-vert.Y(),2) + pow(xyz[2]-vert.Z(),2));
69  if(dist2d>dist) survivingCluster.Add(trackHits[trackHit]);
70  }
71 
72  if(survivingCluster.NCell()>0) gev = survivingCluster.TotalGeV(); // numerator is GeV of all these hits
73 
74  // dx isn't just old track length - dist. To do right would check all places along all trajectory points.
75  // we'll just worry about when one of the trajectory points are too close to vertex. If a track segment went near vertex and it
76  // didn't trigger formation of a new segment (ie a scatter) then it was probably unaffected by vertex activity, so ignore case where
77  // line segment passes near vertex but both trajectory points are outside
78  for(unsigned int p=0;p<track->NTrajectoryPoints()-1;p++) { // to get dx, look through trajectory points and see which are near vertex
79  traj1 = track->TrajectoryPoint(p);
80  traj2 = track->TrajectoryPoint(p+1);
81 
82  if((traj1-vert).Mag() > dist && (traj2-vert).Mag() > dist) length += (traj1-traj2).Mag(); // both end points of segment not near vertex
83  else if((traj1-vert).Mag() < dist && (traj2-vert).Mag() < dist) continue; // both end points too close to vertex, ignore whole segment
84  else { // tricky case of one end point too close and other far away, want to add part of track length that is > dist from vertex
85  TVector3 start, dir; // start at point near vertex, increment away from vertex until we roughly hit the edge of sphere of radius dist
86  if((traj1-vert).Mag() < dist) {
87  start = traj1;
88  dir = (traj2-traj1).Unit();
89  }
90  else {
91  start = traj2;
92  dir = (traj1-traj2).Unit();
93  }
94  TVector3 point = start;
95  bool keepgoing = true;
96  while(keepgoing) { // go down line segment until we hit where we should stop
97  point = point + dir;
98  if( (point - vert).Mag() > (dist - 0.5)) { // this is where we'll stop
99  length += ((traj1-traj2).Mag() - (point-start).Mag()); // add the length left, which is the part farther than dist from vertex
100  keepgoing=false;
101  }
102  if((point-start).Mag() > (traj1-traj2).Mag() + 1) {
103  keepgoing=false;
104  std::cout << "SOMETHING WENT WRONG IN PION RECO DEDX CALC, TOO FAR DISTANCE" << std::endl;
105  }
106  } // end while loop
107  } // end else
108  } // end loop over trajectory points
109  if(length <0) {
110  std::cout << "SOMETHING WRONG IN PION RECO DEDX CALC, LENGTH = " << length << std::endl;
111  return 0;
112  }
113  else if(length == 0) return 0;
114  return gev/length;
115  }
116 
117  //-----------------------------------------------------------------------------------
118  float NDRecoFxs::getTrackActivity(const art::Ptr<rb::Track> track, TVector3 vert, float dist, art::PtrVector< rb::CellHit > const& sliceHits, float activitydist) {
120  TVector3 trackend = track->Stop();
121  // if((track->Start-vert).Mag() > 40 && (track->Stop()-vert).Mag() < 40) trackend = track->Start(); // doesn't do much
122  // in future may want to check start and end positions vs start/end positions of other tracks in slice instead
123 
124  art::PtrVector< rb::CellHit > const& trackHits = track->AllCells();
125  rb::Cluster activityCluster;
126  for (unsigned int sliceHit = 0; sliceHit < sliceHits.size(); ++sliceHit){
127  bool match = 0;
128  for (unsigned int trackHit = 0; trackHit < trackHits.size(); ++trackHit){
129  match = (sliceHits[sliceHit] == trackHits[trackHit]);
130  if (match) break;
131  }//End of loop over hits in track
132 
133  const geo::PlaneGeo* p = geom->Plane(sliceHits[sliceHit]->Plane());
134  double xyz[3] = {0.,0.,0.}, junk = 0;
135  p->Cell(sliceHits[sliceHit]->Cell())->GetCenter(xyz, junk);
136  float disttovert = 0, disttoend = 0;
137  if(sliceHits[sliceHit]->View()==geo::kX) {
138  disttovert = sqrt( pow(xyz[0]-vert.X(),2) + pow(xyz[2]-vert.Z(),2));
139  disttoend = sqrt( pow(xyz[0]-trackend.X(),2) + pow(xyz[2]-trackend.Z(),2));
140  }
141  else if(sliceHits[sliceHit]->View()==geo::kY) {
142  disttovert = sqrt( pow(xyz[1]-vert.Y(),2) + pow(xyz[2]-vert.Z(),2));
143  disttoend = sqrt( pow(xyz[1]-trackend.Y(),2) + pow(xyz[2]-trackend.Z(),2));
144  }
145  if (!match && disttovert>dist && disttoend<activitydist) activityCluster.Add(sliceHits[sliceHit]);
146  }//End of loop over hits in the slice
147  if(activityCluster.NCell()>0) return activityCluster.TotalGeV();
148  else return 0;
149  }
150 
151  //--------------------------------------------------------------------------------------
152 
153  float NDRecoFxs::getTrackProximity(const std::vector< art::Ptr<rb::Track> > sliceTracks, int trknum, TVector3 vert, float dist) {
154 
155  art::Ptr<rb::Track> first = sliceTracks[trknum];
156 
157  TVector3 prox, min = {9999, 9999, 9999}, start1, start2, end1, end2;
158 
159  for(uint i=0; i<sliceTracks.size()-1; i++) {
160 
161  if((int)i==trknum) continue;
162 
163  art::Ptr<rb::Track> second = sliceTracks[i];
164 
165  start1 = first->Start();
166  start2 = second->Start();
167 
168  end1 = first->Stop();
169  end2 = second->Stop();
170 
171  // do all four cases
172  // start with start to start
173 
174  prox = start1-start2;
175  if(prox.Mag() < min.Mag() && (start1-vert).Mag() > dist && (start2-vert).Mag() > dist) min = prox;
176 
177  // start to end
178  prox = start1-end2;
179  if(prox.Mag() < min.Mag() && (start1-vert).Mag() > dist && (end2-vert).Mag() > dist) min = prox;
180 
181  // end to start
182  prox = end1-start2;
183  if(prox.Mag() < min.Mag() && (end1-vert).Mag() > dist && (start2-vert).Mag() > dist) min = prox;
184 
185  // now end to end
186  prox = end1-end2;
187  if(prox.Mag() < min.Mag() && (end1-vert).Mag() > dist && (end2-vert).Mag() > dist) min = prox;
188  }
189  if(min.Mag()>17000) return -1; //no conditions were ever met; min is still {9999,9999,9999)
190  return min.Mag();
191  }
192 
193  //--------------------------------------------------------------------------------------
194  //dedx for prongs is simpler, only one line segment
195  float NDRecoFxs::getProngDedx(const art::Ptr<rb::Prong> prong, TVector3 vert, float dist) {
197  float gev = 0, length = 0;
198  TVector3 traj1, traj2;
199 
200  art::PtrVector< rb::CellHit > const& prongHits = prong->AllCells();
201 
202  rb::Cluster survivingCluster; // make cluster of hits further than dist away from vertex
203 
204  for (unsigned int prongHit = 0; prongHit < prongHits.size(); ++prongHit){ //fill cluster
205  const geo::PlaneGeo* p = geom->Plane(prongHits[prongHit]->Plane());
206  double xyz[3] = {0.,0.,0.}, junk = 0;
207  p->Cell(prongHits[prongHit]->Cell())->GetCenter(xyz, junk);
208  TVector3 cellcenter(xyz);
209  float dist2d = 0;
210  if(prongHits[prongHit]->View()==geo::kX) dist2d = sqrt( pow(xyz[0]-vert.X(),2) + pow(xyz[2]-vert.Z(),2));
211  else if(prongHits[prongHit]->View()==geo::kY) dist2d = sqrt( pow(xyz[1]-vert.Y(),2) + pow(xyz[2]-vert.Z(),2));
212  if(dist2d>dist) survivingCluster.Add(prongHits[prongHit]);
213  }
214 
215  if(survivingCluster.NCell()>0) gev = survivingCluster.TotalGeV(); // numerator is GeV of all these hits
216 
217  // like for tracks but only one line segment
218  traj1 = prong->Start();
219  traj2 = prong->Start() + prong->Dir()*prong->TotalLength();
220 
221  if((traj1-vert).Mag() > dist && (traj2-vert).Mag() > dist) length += (traj1-traj2).Mag(); // both end points of segment not near vertex
222  else if((traj1-vert).Mag() < dist && (traj2-vert).Mag() < dist) length = 0; // both end points too close to vertex, ignore whole segment
223  else { // tricky case of one end point too close and other far away, want to add part of track length that is > dist from vertex
224  TVector3 start, dir; // start at point near vertex, increment away from vertex until we roughly hit the edge of sphere of radius dist
225  if((traj1-vert).Mag() < dist) {
226  start = traj1;
227  dir = (traj2-traj1).Unit();
228  }
229  else {
230  start = traj2;
231  dir = (traj1-traj2).Unit();
232  }
233  TVector3 point = start;
234  bool keepgoing = true;
235  while(keepgoing) { // go down line segment until we hit where we should stop
236  point = point + dir;
237  if( (point - vert).Mag() > (dist - 0.5)) { // this is where we'll stop
238  length += ((traj1-traj2).Mag() - (point-start).Mag()); // add the length left, which is the part farther than dist from vertex
239  keepgoing=false;
240  }
241  if((point-start).Mag() > (traj1-traj2).Mag() + 1) {
242  keepgoing=false;
243  std::cout << "SOMETHING WENT WRONG IN PION RECO DEDX CALC, TOO FAR DISTANCE" << std::endl;
244  }
245  } // end while loop
246  } // end else
247  if(length <0) {
248  std::cout << "SOMETHING WRONG IN PION RECO DEDX CALC, LENGTH = " << length << std::endl;
249  return 0;
250  }
251  else if(length == 0) return 0;
252  return gev/length;
253  }
254 
255  //-----------------------------------------------------------------------------------
256  float NDRecoFxs::getProngActivity(const art::Ptr<rb::Prong> prong, TVector3 vert, float dist, art::PtrVector< rb::CellHit > const& sliceHits, float activitydist) {
258  TVector3 prongend = prong->Start() + prong->Dir()*prong->TotalLength();
259 
260  art::PtrVector< rb::CellHit > const& prongHits = prong->AllCells();
261  rb::Cluster activityCluster;
262  for (unsigned int sliceHit = 0; sliceHit < sliceHits.size(); ++sliceHit){
263  bool match = 0;
264  for (unsigned int prongHit = 0; prongHit < prongHits.size(); ++prongHit){
265  match = (sliceHits[sliceHit] == prongHits[prongHit]);
266  if (match) break;
267  }//End of loop over hits in prong
268 
269  const geo::PlaneGeo* p = geom->Plane(sliceHits[sliceHit]->Plane());
270  double xyz[3] = {0.,0.,0.}, junk = 0;
271  p->Cell(sliceHits[sliceHit]->Cell())->GetCenter(xyz, junk);
272  float disttovert = 0, disttoend = 0;
273  if(sliceHits[sliceHit]->View()==geo::kX) {
274  disttovert = sqrt( pow(xyz[0]-vert.X(),2) + pow(xyz[2]-vert.Z(),2));
275  disttoend = sqrt( pow(xyz[0]-prongend.X(),2) + pow(xyz[2]-prongend.Z(),2));
276  }
277  else if(sliceHits[sliceHit]->View()==geo::kY) {
278  disttovert = sqrt( pow(xyz[1]-vert.Y(),2) + pow(xyz[2]-vert.Z(),2));
279  disttoend = sqrt( pow(xyz[1]-prongend.Y(),2) + pow(xyz[2]-prongend.Z(),2));
280  }
281  if (!match && disttovert>dist && disttoend<activitydist) activityCluster.Add(sliceHits[sliceHit]);
282  }//End of loop over hits in the slice
283  if(activityCluster.NCell()>0) return activityCluster.TotalGeV();
284  else return 0;
285  }
286 
287 
288  //--------------------------------------------------------------------------------------
289  float NDRecoFxs::getProngProximity(const std::vector< art::Ptr<rb::Prong> > sliceProngs, int pngnum, TVector3 vert, float dist) {
290 
291  art::Ptr<rb::Prong> first = sliceProngs[pngnum];
292 
293  TVector3 prox, min = {9999, 9999, 9999}, start1, start2, end1, end2;
294 
295  for(uint i=0; i<sliceProngs.size()-1; i++) {
296 
297  if((int)i==pngnum) continue;
298 
299  art::Ptr<rb::Prong> second = sliceProngs[i];
300 
301  start1 = first->Start();
302  start2 = second->Start();
303 
304  end1 = first->Start() + first->Dir()*first->TotalLength();
305  end2 = second->Start() + first->Dir()*first->TotalLength();
306 
307  // do all four cases
308  // start with start to start
309 
310  prox = start1-start2;
311  if(prox.Mag() < min.Mag() && (start1-vert).Mag() > dist && (start2-vert).Mag() > dist) min = prox;
312 
313  // start to end
314  prox = start1-end2;
315  if(prox.Mag() < min.Mag() && (start1-vert).Mag() > dist && (end2-vert).Mag() > dist) min = prox;
316 
317  // end to start
318  prox = end1-start2;
319  if(prox.Mag() < min.Mag() && (end1-vert).Mag() > dist && (start2-vert).Mag() > dist) min = prox;
320 
321  // now end to end
322  prox = end1-end2;
323  if(prox.Mag() < min.Mag() && (end1-vert).Mag() > dist && (end2-vert).Mag() > dist) min = prox;
324  }
325  if(min.Mag()>17000) return -1; //no conditions were ever met; min is still {9999,9999,9999)
326  return min.Mag();
327  }
328 
329 
330 } //end namespace
size_t NTrajectoryPoints() const
Definition: Track.h:83
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
float getProngActivity(const art::Ptr< rb::Prong > prong, TVector3 vert, float dist, art::PtrVector< rb::CellHit > const &sliceHits, float activitydist)
Definition: NDRecoFxs.cxx:256
float getTrackProximity(const std::vector< art::Ptr< rb::Track > > sliceTracks, int trknum, TVector3 vertex, float dist)
Definition: NDRecoFxs.cxx:153
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
float getProngProximity(const std::vector< art::Ptr< rb::Prong > > sliceProngs, int pngnum, TVector3 vertex, float dist)
Definition: NDRecoFxs.cxx:289
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
float getTrackScatt(const art::Ptr< rb::Track > track, TVector3 vert, float dist)
Definition: NDRecoFxs.cxx:32
float getTrackDedx(const art::Ptr< rb::Track > track, TVector3 vert, float dist)
Definition: NDRecoFxs.cxx:53
T acos(T number)
Definition: d0nt_math.hpp:54
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
double dist
Definition: runWimpSim.h:113
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Definition: Prong.cxx:186
length
Definition: demo0.py:21
olddir
Are we streaming the files via xrootd? txtcmd="cat %s | xargs -n1 samweb2xrootd > xrootd_inFile...
Definition: runNovaSAM.py:800
float getProngDedx(const art::Ptr< rb::Prong > prong, TVector3 vert, float dist)
Definition: NDRecoFxs.cxx:195
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
Collect Geo headers and supply basic geometry functions.
TVector3 Unit() const
Definition: View.py:1
size_type size() const
Definition: PtrVector.h:308
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
float getTrackActivity(const art::Ptr< rb::Track > track, TVector3 vert, float dist, art::PtrVector< rb::CellHit > const &sliceHits, float activitydist)
Definition: NDRecoFxs.cxx:118
TDirectory * dir
Definition: macro.C:5
void geom(int which=0)
Definition: geom.C:163
float Mag() const
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t e
Definition: plot.C:35
produce variables useful to reconstruction pion energy
Definition: FillReco.h:12
virtual ~NDRecoFxs()
Definition: NDRecoFxs.cxx:26
Encapsulate the geometry of one entire detector (near, far, ndos)
unsigned int uint