TrackInfoFxs.cxx
Go to the documentation of this file.
2 
4 #include "Geometry/Geometry.h"
6 #include "GeometryObjects/Geo.h"
8 #include "RecoBase/RecoHit.h"
9 
11 
12 #include "TMath.h"
13 #include "TVector3.h"
14 
15 #include "MCCheater/BackTracker.h"
17 
18 namespace trackinfo
19 {
20 
21  //-------------------------------------------------------
23  }
24 
25  //-------------------------------------------------------
27  }
28 
29  //-----------------------------------------------------------------------------------
30 
31  float TrackInfoFxs::getAveTrackdEdxLastCells( const art::Ptr<rb::Track> track, unsigned int nLastCells){
32  // This function returns the average dEdx in MeV/cm for the last nLastCells cells
33 
34  float avededxlast = -5; //default value
35 
36  if(track->NCell() < nLastCells){ //the number of cells in the track can not be less than the required number of last cells
37  return avededxlast;
38  }
39  else{
41 
42  std::vector<std::vector<double> > vec_trkRecoHitInfo;
43  std::pair<uint32_t, uint32_t> pc;
44  TVector3 point;
45  std::set<double> planeZBounds;
46  calib::FindZBoundaries(planeZBounds);
47  auto boundaryMap = calib::MakeZBoundaryMap(planeZBounds, track->Trajectory());
48  if(boundaryMap.size() < 1) return avededxlast;
49  //Loop over muon track hits
50  for(unsigned int jcell = 0; jcell < track->NCell();jcell++){
51  const art::Ptr<rb::CellHit> cellhit = track->Cell(jcell);
52  if(btrk->IsNoise(cellhit)) continue;
53  pc.first = cellhit->Plane();
54  pc.second = cellhit->Cell();
55 
56  rb::RecoHit rhit = track->RecoHit(cellhit);
57  if(!rhit.IsCalibrated()) continue;
58 
59  std::vector<double> rhitInfo(5);
60  rhitInfo[0] = rhit.X();
61  rhitInfo[1] = rhit.Y();
62  rhitInfo[2] = rhit.Z();
63  point = TVector3(rhitInfo[0], rhitInfo[1], rhitInfo[2]);
64  rhitInfo[3] = rhit.GeV();
65  rhitInfo[4] = calib::PathLengthInCell(boundaryMap,point,pc);
66 
67  vec_trkRecoHitInfo.push_back(rhitInfo);
68  } //end loop over trk cell hits
69 
70  // using lambda function to sort it
71  std::sort(vec_trkRecoHitInfo.begin(), vec_trkRecoHitInfo.end(),
72  [](const std::vector<double>& a, const std::vector<double>& b) {
73  return a[2] < b[2];
74  });
75 
76  float sumEnLast = 0.0;
77  float sumPathLast = 0.0;
78  if (!vec_trkRecoHitInfo.empty()){
79  if (nLastCells < 1 || vec_trkRecoHitInfo.size() < nLastCells){
80  mf::LogWarning("TrackInfoFxs") << "nLastCells incompatible with reco hit info size";
81  return avededxlast;
82  }
83 
84  size_t min_index = vec_trkRecoHitInfo.size()-nLastCells;
85  for(size_t a = min_index;a<vec_trkRecoHitInfo.size();a++) {
86  sumEnLast += (1000.*vec_trkRecoHitInfo[a][3]); // From GeV to MeV
87  sumPathLast += vec_trkRecoHitInfo[a][4]; // Distance in cm
88  }
89  }
90 
91  if(sumPathLast > 0){
92  avededxlast = (sumEnLast)/(sumPathLast); // average dEdx in MeV/cm
93  }
94  }
95  return avededxlast;
96  }//end of getAveTrackdEdxLastCells
97 
98  float TrackInfoFxs::getAveTrackdEdxEnd(std::vector<double> &dEtemp, std::vector<double> &dxtemp,
99  std::vector<double> &stemp, double trk_len, double wanted_dist){
100  float avededxend = -5; //default value
101  if(trk_len < wanted_dist){
102  return avededxend;
103  }
104 
105  double acc_dE = 0;
106  double acc_dx = 0;
107  //finding the minimum index:
108  unsigned int min_index = 10000;
109  for(unsigned int p = 0; p< stemp.size(); p++){
110  double distcell_to_end = trk_len-stemp[p];
111  if(distcell_to_end < 0){
112  mf::LogWarning("TrackInfoFxs")<<"Something is wrong with this track, distance of a cell in the track to the end of the track is negative";
113  return avededxend;
114  }
115  if(distcell_to_end < wanted_dist){
116  min_index = p;
117  break;
118  }
119  }
120 
121  for(unsigned int p = min_index; p < stemp.size(); p++){
122  if(dxtemp[p] <= 0){
123  mf::LogWarning("TrackInfoFxs")<<"Something is wrong with dx... not positive value";
124  return avededxend;
125  }
126  acc_dE += 1000. * dEtemp[p]; //in MeV
127  acc_dx += dxtemp[p]; //in cm
128  }
129  if(acc_dx>0)avededxend = acc_dE / acc_dx;
130 
131  return avededxend;
132  }//end of getAveTrackdEdxEnd
133 
134 } //end namespace
back track the reconstruction to the simulation
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
zBoundMap MakeZBoundaryMap(std::set< double > const &planeZBounds, std::vector< TVector3 > const &trajectory)
Return a map of the z position of each trajectory point on a track to the bounding positions of where...
Definition: CalibUtil.cxx:354
unsigned short Plane() const
Definition: CellHit.h:39
const char * p
Definition: xmltok.h:285
float Z() const
Definition: RecoHit.h:38
Definition: event.h:19
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
unsigned short Cell() const
Definition: CellHit.h:40
const double a
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
Collect Geo headers and supply basic geometry functions.
float getAveTrackdEdxEnd(std::vector< double > &dEtemp, std::vector< double > &dxtemp, std::vector< double > &stemp, double trk_len, double wanted_dist)
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
float getAveTrackdEdxLastCells(const art::Ptr< rb::Track > track, unsigned int nLastCells)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
Track information.
Definition: FillReco.h:15
float GeV() const
Definition: RecoHit.cxx:69
float X() const
Definition: RecoHit.h:36
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
Definition: Track.h:86
const hit & b
Definition: hits.cxx:21
float Y() const
Definition: RecoHit.h:37
void FindZBoundaries(std::set< double > &planeZBounds)
Find the boundaries in the z direction of planes in the detector.
Definition: CalibUtil.cxx:332
Encapsulate the geometry of one entire detector (near, far, ndos)
float PathLengthInCell(zBoundMap const &zBounds, TVector3 const &recoHitLoc, std::pair< uint32_t, uint32_t > const &pc)
Return the path length of a track in the cell in question.
Definition: CalibUtil.cxx:498