RegPixelMapProducer.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file RegPixelMapProducer.h
3 /// \brief RegPixelMapProducer for CVN
4 /// \author Dominick Rocco - rocco@physics.umn.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include <iostream>
8 #include <ostream>
9 #include <list>
10 #include <algorithm>
11 
13 #include "CVN/func/AssignLabels.h"
14 #include "TVector2.h"
15 #include "RecoBase/CellHit.h"
16 #include "RecoBase/RecoHit.h"
17 
18 
19 namespace cvn
20 {
21 
22  RegPixelMapProducer::RegPixelMapProducer(unsigned int nPlane, unsigned int nCell,
23  unsigned int maxPlaneGap, bool useGeV):
24  fNPlane(nPlane),
25  fNCell(nCell),
26  fMaxPlaneGap(maxPlaneGap),
27  fUseGeV(useGeV)
28  {}
29 
30 
32  {
33 
34  Boundary bound = DefineBoundary(cluster);
35 
36  return CreateMapGivenBoundary(cluster, bound);
37 
38 
39 
40  }
41 
43  const Boundary& bound)
44  {
45 
46  RegPixelMap pm(fNPlane, fNCell, bound);
47  for(size_t iHit = 0; iHit < cluster.NCell(); ++iHit)
48  {
49  unsigned int plane = cluster.Cell(iHit)->Plane();
50  unsigned int cell = cluster.Cell(iHit)->Cell();
51  rb::RecoHit recoHit = cluster.RecoHit(iHit);
52 
53  art::Ptr< rb::CellHit >cellHit = cluster.Cell(iHit);
54 
55  float value;
56  if (fUseGeV) {
57  // Use GeV if we've got it. Otherwise use GeV equivalent of PE.
58  float gev;
59  if(recoHit.IsCalibrated())
60  gev = recoHit.GeV();
61  else
62  { // Put something rough in there, just as a placeholder.
63  gev = cluster.Cell(iHit)->PE() / 15123.7; // Rough PE->GeV scale.
64  }
65 
66  value = gev;
67  }
68  else {
69  // Use PECorr if we've got it. Otherwise use PE.
70  float pe = recoHit.IsCalibrated()? recoHit.PECorr() :
71  cluster.Cell(iHit)->PE();
72  value = pe;
73  }
74 
75  double purity =0.0;
77  label = HitClassify( cellHit, &label, &purity);
78  pm.Add(plane, cell, value, label, purity);
79  }
80 // std::cout<<"nHit "<<cluster.NCell()<<std::endl;
81  return pm;
82 
83  }
84 
85  RegPixelMap RegPixelMapProducer::CreateMapGivenShowerVertex(const rb::Shower& shower, const Boundary& bound, unsigned int vtxPlane, unsigned int vtxXCell, unsigned int vtxYCell)
86  {
87  RegPixelMap pm(fNPlane, fNCell, vtxPlane, vtxXCell, vtxYCell, bound);
88  for(size_t iHit = 0; iHit < shower.NCell(); ++iHit)
89  {
90  const unsigned int plane = shower.Cell(iHit)->Plane();
91  const unsigned int cell = shower.Cell(iHit)->Cell();
92  rb::RecoHit recoHit = shower.RecoHit(iHit);
93  geo::View_t view = shower.Cell(iHit)->View();
94  int cellView = -1;
95  if(view == geo::kX){
96  cellView = 0;
97  }else if(view == geo::kY){
98  cellView = 1;
99  }
100 
101  art::Ptr< rb::CellHit >cellHit = shower.Cell(iHit);
102 
103 
104  float gev;
105  if(recoHit.IsCalibrated())
106  gev = recoHit.GeV();
107  else
108  {
109  gev = shower.Cell(iHit)->PE() / 15123.7;
110  }
111 // const float value = shower.HitDeconvE(iHit);
112 // const float value = gev*shower.Weight(iHit);;
113  const float value = gev;
114  double purity =0.0;
116  label = HitClassify( cellHit, &label, &purity);
117 // pm.Add(plane, cell, value, label, purity);
118 //std::cout<<"shower.size, plane, cell, cellView, value, fNPlane, fNCell "<<shower.NCell()<<" "<<plane<<" "<<cell<<" "<<cellView<<" "<<value<<" "<<fNPlane<<" "<<fNCell<<std::endl;
119 // pm.Add(plane, cell, cellView, value, label, purity, 100.);//prong
120  pm.Add(plane, cell, cellView, value, label, purity);//prong
121  }
122 // std::cout<<"nHit "<<shower.NCell()<<std::endl;
123  return pm;
124 
125  }
126 
127 
128  RegPixelMap RegPixelMapProducer::CreateMapGivenVertex(const rb::Cluster& cluster, const Boundary& bound, unsigned int vtxPlane, unsigned int vtxXCell, unsigned int vtxYCell){
129  RegPixelMap pm(fNPlane, fNCell, vtxPlane, vtxXCell, vtxYCell, bound);
130 // int nhit=0;
131 // float cle=0.;
132 
133  for(size_t iHit = 0; iHit < cluster.NCell(); ++iHit)
134  {
135  unsigned int plane = cluster.Cell(iHit)->Plane();
136  unsigned int cell = cluster.Cell(iHit)->Cell();
137  rb::RecoHit recoHit = cluster.RecoHit(iHit);
138  geo::View_t view = cluster.Cell(iHit)->View();
139  int cellView = -1;
140  if(view == geo::kX){
141  cellView = 0;
142  }else if(view == geo::kY){
143  cellView = 1;
144  }
145  art::Ptr< rb::CellHit >cellHit = cluster.Cell(iHit);
146 
147  float value;
148  if (fUseGeV) {
149  // Use GeV if we've got it. Otherwise use GeV equivalent of PE.
150  float gev;
151  if(recoHit.IsCalibrated())
152  gev = recoHit.GeV();
153  else
154  { // Put something rough in there, just as a placeholder.
155  gev = cluster.Cell(iHit)->PE() / 15123.7; // Rough PE->GeV scale.
156  }
157 
158  value = gev;
159  }
160  else {
161  // Use PECorr if we've got it. Otherwise use PE.
162  float pe = recoHit.IsCalibrated()? recoHit.PECorr() :
163  cluster.Cell(iHit)->PE();
164  value = pe;
165  }
166 
167  double purity =0.0;
169  label = HitClassify( cellHit, &label, &purity);
170 // pm.Add(plane, cell, value, label, purity);
171 //std::cout<<"cluster.size, plane, cell, cellView, value, fNPlane, fNCell "<<cluster.NCell()<<" "<<plane<<" "<<cell<<" "<<cellView<<" "<<value<<" "<<fNPlane<<" "<<fNCell<<std::endl;
172  pm.Add(plane, cell, cellView, value, label, purity);
173 // nhit++;
174  }
175 // std::cout<<"Cluster Event nHit cluster "<<nhit<<std::endl;
176 // std::cout<<"Cluster Event Energy cle "<<" "<<cle<<std::endl;
177  return pm;
178  }
179 
180 
181 
182  RegPixelMap RegPixelMapProducer::CreateMapGivenVertex(const std::vector<art::Ptr<rb::Shower>> showers3D, const Boundary& bound, unsigned int vtxPlane, unsigned int vtxXCell, unsigned int vtxYCell)
183  {
184 
185  RegPixelMap pm(fNPlane, fNCell, vtxPlane, vtxXCell, vtxYCell, bound);
186 // int nhit=0;
187  for( unsigned int iProng = 0; iProng < showers3D.size(); ++iProng ){
188  for(size_t iHit = 0; iHit < showers3D[iProng]->NCell(); ++iHit)
189  {
190  const unsigned int plane = showers3D[iProng]->Cell(iHit)->Plane();
191  const unsigned int cell = showers3D[iProng]->Cell(iHit)->Cell();
192  rb::RecoHit recoHit = showers3D[iProng]->RecoHit(iHit);
193  geo::View_t view = showers3D[iProng]->Cell(iHit)->View();
194  int cellView = -1;
195  if(view == geo::kX){
196  cellView = 0;
197  }else if(view == geo::kY){
198  cellView = 1;
199  }
200 
201  art::Ptr< rb::CellHit >cellHit = showers3D[iProng]->Cell(iHit);
202 
203  float gev;
204  if(recoHit.IsCalibrated())
205  gev = recoHit.GeV();
206  else
207  {
208  gev = showers3D[iProng]->Cell(iHit)->PE() / 15123.7;
209  }
210 // const float value = showers3D[iProng]->HitDeconvE(iHit);
211  const float value = gev;
212 
213  double purity =0.0;
215  label = HitClassify( cellHit, &label, &purity);
216  pm.Add(plane, cell, cellView, value, label, purity);
217 // nhit++;
218  }
219  }
220 // std::cout<<"Shower Event nHit sh "<<nhit<<std::endl;
221  return pm;
222  }
223 
224 
225 
226  std::ostream& operator<<(std::ostream& os, const RegPixelMapProducer& p)
227  {
228  os << "RegPixelMapProducer: "
229  << p.NCell() <<" cells X " << p.NPlane() << " planes";
230  return os;
231  }
232 
233 
235  std::set<unsigned int> planes;
236  for(size_t iHit = 0; iHit < cluster.NCell(); ++iHit)
237  planes.insert(cluster.Cell(iHit)->Plane());
238 
239  unsigned int firstPlane = cluster.MinPlane();
240  unsigned int previousPlane = cluster.MinPlane();
241  for(const auto& plane:planes){
242  unsigned int gap = plane - previousPlane;
243  if (gap > fMaxPlaneGap)
244  {
245  firstPlane = plane;
246  }
247  previousPlane = plane;
248  }
249  return firstPlane;
250  }
251 
253  const rb::Cluster& cluster)
254  {
255 
256  unsigned int window = 8;
257  unsigned int threshold = 4;
258 
259  std::map<unsigned int, unsigned int> planeMap;
260 
261  /// Count the hits in each plane
262  for(size_t iHit = 0; iHit < cluster.NCell(); ++iHit)
263  planeMap[cluster.Cell(iHit)->Plane()]++;
264 
265 
266 
267 
268  for(unsigned int iPlane = cluster.MinPlane();
269  iPlane < cluster.MaxPlane() - window; ++iPlane)
270  {
271  // Don't start a window unless there's a hit in this plane
272  if(not planeMap.count(iPlane)) continue;
273 
274  // Loop over planes in the window to count up hits
275  unsigned int windowCount = 0;
276  for(unsigned int iWindow = iPlane; iWindow - iPlane < window; ++iWindow)
277  {
278  /// If we have hits in that plane, count them.
279  if(planeMap.count(iWindow))
280  windowCount += planeMap[iWindow];
281  }
282 
283  // If our window count is above threshold, this is the plane we want
284  if(windowCount >= threshold)
285  return iPlane;
286 
287  }
288  // If we're here, we didn't find it. Punt by returning min plane.
289  return cluster.MinPlane();
290 
291  }
292 
293 
294  std::array<unsigned int, 2> RegPixelMapProducer::FindCenterMedian(
295  const rb::Cluster& cluster,
296  unsigned int firstPlane)
297  {
298 
299  std::vector<unsigned int> cells[2];
300  for(unsigned int iHit = 0; iHit < cluster.NCell(); ++iHit)
301  {
302  unsigned int plane = cluster.Cell(iHit)->Plane();
303  unsigned int view = plane%2;
304  if(not InPlaneRcvne(plane, firstPlane)) continue;
305  /// add this cell to the list for the appropriate view
306  cells[view].push_back(cluster.Cell(iHit)->Cell());
307  }
308  // Containers for and modes and maxima as we loop
309  std::array<unsigned int, 2> median;
310 
311  // Do views independently
312  for(unsigned int view = 0; view < 2; ++view)
313  {
314  std::sort(cells[view].begin(), cells[view].end());
315  if(cells[view].size() > 0)
316  median[view] = cells[view][cells[view].size()/2];
317  else
318  median[view] = 0 ; // It's ok to fall back since there are no hits
319 
320  }
321  return median;
322 
323 
324 
325  }
326 
328  {
329 
330  unsigned int firstPlane = FindVertexWindowThreshold(cluster);
331  std::array<unsigned int, 2> centerCell = FindCenterMedian(cluster,
332  firstPlane);
333 
334  Boundary bound(fNPlane, fNCell, firstPlane,
335  centerCell[0], centerCell[1]);
336 
337 
338  return bound;
339  }
340 
341 
342 
344  unsigned int firstPlane)
345  {
346  return plane >= firstPlane && plane < firstPlane + fNPlane;
347  }
348 
349  bool RegPixelMapProducer::InCellRcvne(unsigned int cell, float meanCell){
350  return fabs(cell - meanCell) <= (float) (fNCell / 2);
351  }
352 
353 
354 }
unsigned int NCell() const
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
Boundary DefineBoundary(const rb::Cluster &cluster)
Get boundaries for pixel map representation of cluster.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Producer algorithm for RegPixelMap, input to CVN neural net.
enum cvn::HType HitType
void Add(const unsigned int &plane, const unsigned int &cell, const float &pe, const HitType label, const double purity)
Definition: RegPixelMap.cxx:65
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::ostream & operator<<(std::ostream &os, const PixelMapProducer &p)
A collection of associated CellHits.
Definition: Cluster.h:47
Defines an enumeration for prong classification.
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
const char * label
unsigned short Cell() const
Definition: CellHit.h:40
RegPixelMap CreateMapGivenBoundary(const rb::Cluster &cluster, const Boundary &bound)
RegPixelMap CreateMapGivenShowerVertex(const rb::Shower &shower, const Boundary &bound, unsigned int vtxPlane, unsigned int vtxXCell, unsigned int vtxYCell)
unsigned int fNCell
Number of cells, width of pixel map.
const XML_Char int const XML_Char * value
Definition: expat.h:331
bool InCellRcvne(unsigned int cell, float meanCell)
Determine whether cell is within cell rcvne based on meanCell.
RegPixelMapProducer for CVN.
unsigned int FindVertexWindowThreshold(const rb::Cluster &cluster)
Find vertex plane by requiring gaps are not wider than maxGap.
float PE() const
Definition: CellHit.h:42
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
std::array< unsigned int, 2 > FindCenterMedian(const rb::Cluster &cluster, unsigned int firstPlane)
Find center cell based on mean of hits in each view (2 views)
double median(TH1D *h)
Definition: absCal.cxx:524
RegPixelMap, basic input to CVN neural net.
Definition: RegPixelMap.h:22
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
RegPixelMap CreateMapGivenVertex(const rb::Cluster &cluster, const Boundary &bound, unsigned int vtxPlane, unsigned int vtxXCell, unsigned int vtxYCell)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
RegPixelMapProducer(unsigned int nPlane, unsigned int nCell, unsigned int maxPlaneGap, bool useGeV)
float GeV() const
Definition: RecoHit.cxx:69
A rb::Prong with a length.
Definition: Shower.h:18
unsigned int fNPlane
Number of planes, length for pixel maps.
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
RegPixelMap CreateMap(const rb::Cluster &slice)
float PECorr() const
Definition: RecoHit.cxx:47
HitType HitClassify(art::Ptr< rb::CellHit > Hit, HitType *hType, double *hPurity)
bool InPlaneRcvne(unsigned int plane, unsigned int firstPlane)
Determine plane is within plane rcvne of a PixelMap.
unsigned int NPlane() const
unsigned int FindVertexMaxGap(const rb::Cluster &cluster)
Find vertex plane by requiring gaps are not wider than maxGap.