PixelMapProducer.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file PixelMapProducer.h
3 /// \brief PixelMapProducer 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  PixelMapProducer::PixelMapProducer(unsigned int nPlane, unsigned int nCell,
23  unsigned int maxPlaneGap, bool useGeV,
24  bool useFullND, bool hitsb, bool timeb,
25  bool fullSlice, bool fillHitNu,
26  bool PEOnly):
27  fNPlane(nPlane),
28  fNCell(nCell),
29  fMaxPlaneGap(maxPlaneGap),
30  fUseGeV(useGeV),
31  fUseFullND(useFullND),
32  fHits(hitsb),
33  fTime(timeb),
34  fFullSlice(fullSlice),
35  fFillHitNu(fillHitNu),
36  fPEOnly(PEOnly)
37  {}
38 
39 
40  ////////////////////////////////////////////////////////////////////////
42  {
43 
44  Boundary bound = DefineBoundary(cluster);
45 
46  return CreateMapGivenBoundary(cluster, bound);
47 
48 
49 
50  }
51 
52  ////////////////////////////////////////////////////////////////////////
54  const Boundary& bound)
55  {
56 
57  PixelMap pm(fNPlane, fNCell, bound, fPEOnly);
58 
59  std::map<int,int> objectMap;
60  int objectCount = 0;
61 
62  for(size_t iHit = 0; iHit < cluster.NCell(); ++iHit)
63  {
64  art::Ptr< rb::CellHit >cellHit = cluster.Cell(iHit);
65  rb::RecoHit recoHit = cluster.RecoHit(iHit);
66  geo::View_t view = cluster.Cell(iHit)->View();
67 
68  // Get values to fill PixelMap with
69  unsigned int view_uint = (unsigned int)view;
70  unsigned int hit_uint = (unsigned int)iHit;
71  unsigned int plane = cluster.Cell(iHit)->Plane();
72  unsigned int cell = cluster.Cell(iHit)->Cell();
73  float hitPE = cluster.Cell(iHit)->PE();
74  float time = cluster.Cell(iHit)->TNS();
75  float corr = 15123.7; // Rough PE to GeV? Rev 17083
76  float gev = ( recoHit.IsCalibrated() ) ? recoHit.GeV() : hitPE/corr;
77  float pe = ( recoHit.IsCalibrated() ) ? recoHit.PECorr() : hitPE;
78  float value = ( fUseGeV ) ? gev : pe;
79  double purity = 0.0;
81  int object = 0;
82  int neutrino = 0;
83 
84  // Get Hit ID
85  label = HitClassify( cellHit, &label, &purity);
86  // Get Object ID
87  object = HitObject( cellHit, objectCount, objectMap);
88  // Get Neutrino Index
89  if (fFillHitNu) {
90  neutrino = HitNuIndex( cellHit, cluster );
91  }
92  else{
93  neutrino = -5;
94  }
95 
96  // Overwrite PixelMap main value for time or hit? (this is
97  // from revision 26153, probably should find a way to avoid
98  // fundamentally changing what an ART object like PixelMaps
99  // contain without warning to the user of the ART file.
100  if (fHits) value = 1.0;
101  else if (fTime) value = time;
102 
103  pm.Add(plane, cell, value, view_uint, hit_uint, time, label, purity, object, neutrino);
104  }
105 
106  pm.fFracX = (float) pm.fOccupiedPixelsX / (float) cluster.NXCell();
107  pm.fFracY = (float) pm.fOccupiedPixelsY / (float) cluster.NYCell();
108 
109  return pm;
110 
111  }
112 
113 
114  ////////////////////////////////////////////////////////////////////////
115  std::ostream& operator<<(std::ostream& os, const PixelMapProducer& p)
116  {
117  os << "PixelMapProducer: "
118  << p.NCell() <<" cells X " << p.NPlane() << " planes";
119  return os;
120  }
121 
122 
123  ////////////////////////////////////////////////////////////////////////
124  unsigned int PixelMapProducer::FindVertexMaxGap(const rb::Cluster& cluster){
125  std::set<unsigned int> planes;
126  for(size_t iHit = 0; iHit < cluster.NCell(); ++iHit)
127  planes.insert(cluster.Cell(iHit)->Plane());
128 
129  unsigned int firstPlane = cluster.MinPlane();
130  unsigned int previousPlane = cluster.MinPlane();
131  for(const auto& plane:planes){
132  unsigned int gap = plane - previousPlane;
133  if (gap > fMaxPlaneGap)
134  {
135  firstPlane = plane;
136  }
137  previousPlane = plane;
138  }
139  return firstPlane;
140  }
141 
142  ////////////////////////////////////////////////////////////////////////
144  const rb::Cluster& cluster)
145  {
146 
147  unsigned int window = 8;
148  unsigned int threshold = 4;
149 
150  std::map<unsigned int, unsigned int> planeMap;
151 
152  /// Count the hits in each plane
153  for(size_t iHit = 0; iHit < cluster.NCell(); ++iHit)
154  planeMap[cluster.Cell(iHit)->Plane()]++;
155 
156 
157  for(unsigned int iPlane = cluster.MinPlane();
158  iPlane < cluster.MaxPlane() - window; ++iPlane)
159  {
160  // Don't start a window unless there's a hit in this plane
161  if(not planeMap.count(iPlane)) continue;
162 
163  // Loop over planes in the window to count up hits
164  unsigned int windowCount = 0;
165  for(unsigned int iWindow = iPlane; iWindow - iPlane < window; ++iWindow)
166  {
167  /// If we have hits in that plane, count them.
168  if(planeMap.count(iWindow))
169  windowCount += planeMap[iWindow];
170  }
171 
172  // If our window count is above threshold, this is the plane we want
173  if(windowCount >= threshold)
174  return iPlane;
175 
176  }
177  // If we're here, we didn't find it. Punt by returning min plane.
178  return cluster.MinPlane();
179 
180  }
181 
182 
183  ////////////////////////////////////////////////////////////////////////
184  std::array<unsigned int, 2> PixelMapProducer::FindCenterMedian(
185  const rb::Cluster& cluster,
186  unsigned int firstPlane)
187  {
188 
189  std::vector<unsigned int> cells[2];
190  for(unsigned int iHit = 0; iHit < cluster.NCell(); ++iHit)
191  {
192  unsigned int plane = cluster.Cell(iHit)->Plane();
193  geo::View_t view = cluster.Cell(iHit)->View();
194 
195  if(not InPlaneRcvne(plane, firstPlane)) continue;
196  /// add this cell to the list for the appropriate view
197  cells[view].push_back(cluster.Cell(iHit)->Cell());
198  }
199  // Containers for and modes and maxima as we loop
200  std::array<unsigned int, 2> median;
201 
202  // Do views independently
203  for(unsigned int view = 0; view < 2; ++view)
204  {
205  std::sort(cells[view].begin(), cells[view].end());
206  if(cells[view].size() > 0)
207  median[view] = cells[view][cells[view].size()/2];
208  else
209  median[view] = 0 ; // It's ok to fall back since there are no hits
210 
211  }
212  return median;
213 
214  }
215 
216  ////////////////////////////////////////////////////////////////////////
218  {
219  unsigned int firstPlane;
220  unsigned int firstView;
221  std::array<unsigned int, 2> centerCell;
222  if (fFullSlice){
223  fNPlane = cluster.MaxPlane() - cluster.MinPlane() + 1;
224  fNPlane += fNPlane % 2;
225  firstPlane = cluster.MinPlane();
226  firstView = (firstPlane +1) % 2;
227  fNCell = cluster.MaxCell(geo::kXorY) - cluster.MinCell(geo::kXorY) + 1;
228  fNCell += fNCell % 2;
229  centerCell[0] = cluster.MinCell(geo::kXorY) + fNCell/2;
230  centerCell[1] = cluster.MinCell(geo::kXorY) + fNCell/2;
231  }
232  else{
233  firstPlane = FindVertexWindowThreshold(cluster);
234  //TODO: Need to store view of first plane for pixel map drawing later. In ND and FD Y view is even planes, X view is odd. X = view 0, Y = view 1. Better to use the geometry service, but debating if it is worth the dependency for one line. There may be a better way. Putting a note here and will come back
235  firstView = (firstPlane +1) % 2;
236  if(fUseFullND)
237  firstPlane = 0;
238  centerCell = FindCenterMedian(cluster,firstPlane);
239  }
240 
241  Boundary bound(fNPlane, fNCell, firstPlane,
242  centerCell[0], centerCell[1], firstView);
243 
244 
245  return bound;
246  }
247 
248 
249 
250  ////////////////////////////////////////////////////////////////////////
252  unsigned int firstPlane)
253  {
254  return plane >= firstPlane && plane < firstPlane + fNPlane;
255  }
256 
257  ////////////////////////////////////////////////////////////////////////
258  bool PixelMapProducer::InCellRcvne(unsigned int cell, float meanCell){
259  return fabs(cell - meanCell) <= (float) (fNCell / 2);
260  }
261 
262 
263 }
float TNS() const
Definition: CellHit.h:46
unsigned int fOccupiedPixelsX
Number of non-zero entries in pixel map. Counted when the Add function is called. ...
Definition: PixelMap.h:124
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
int HitObject(art::Ptr< rb::CellHit > Hit, int &hObjectCount, std::map< int, int > &hObjectMap)
enum cvn::HType HitType
X or Y views.
Definition: PlaneGeo.h:30
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
int HitNuIndex(art::Ptr< rb::CellHit > Hit, const rb::Cluster &evtCluster)
std::ostream & operator<<(std::ostream &os, const PixelMapProducer &p)
A collection of associated CellHits.
Definition: Cluster.h:47
unsigned int MaxCell(geo::View_t view) const
Definition: Cluster.cxx:518
double corr(TGraph *g, double thresh)
void Add(const unsigned int &plane, const unsigned int &cell, const float &pe, unsigned int &view, unsigned int &hitID, const float &time, const HitType label, const double purity, const int object, const int nuIdx)
Definition: PixelMap.cxx:65
Defines an enumeration for prong classification.
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
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
PixelMapProducer(unsigned int nPlane, unsigned int nCell, unsigned int maxPlaneGap, bool useGeV, bool useFullND, bool hits, bool time, bool fullSlice=false, bool fillHitNu=false, bool PEOnly=false)
unsigned int NCell() const
const XML_Char int const XML_Char * value
Definition: expat.h:331
Boundary DefineBoundary(const rb::Cluster &cluster)
Get boundaries for pixel map representation of cluster.
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
unsigned int fOccupiedPixelsY
Number of non-zero entries in pixel map. Counted when the Add function is called. ...
Definition: PixelMap.h:125
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
float fFracY
fraction of y view hits contained in map
Definition: PixelMap.h:120
unsigned int FindVertexMaxGap(const rb::Cluster &cluster)
Find vertex plane by requiring gaps are not wider than maxGap.
PixelMap CreateMapGivenBoundary(const rb::Cluster &cluster, const Boundary &bound)
double median(TH1D *h)
Definition: absCal.cxx:524
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)
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
unsigned int fNPlane
Number of planes, length for pixel maps.
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
bool InPlaneRcvne(unsigned int plane, unsigned int firstPlane)
Determine plane is within plane rcvne of a PixelMap.
float GeV() const
Definition: RecoHit.cxx:69
float fFracX
fraction of x view hits contained in map
Definition: PixelMap.h:119
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
unsigned int MinCell(geo::View_t view) const
Definition: Cluster.cxx:472
PixelMap, basic input to CVN neural net.
Definition: PixelMap.h:23
PixelMapProducer for CVN.
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
bool InCellRcvne(unsigned int cell, float meanCell)
Determine whether cell is within cell rcvne based on meanCell.
float PECorr() const
Definition: RecoHit.cxx:47
HitType HitClassify(art::Ptr< rb::CellHit > Hit, HitType *hType, double *hPurity)
unsigned int NPlane() const
PixelMap CreateMap(const rb::Cluster &slice)
Producer algorithm for PixelMap, input to CVN neural net.
unsigned int FindVertexWindowThreshold(const rb::Cluster &cluster)
Find vertex plane by requiring gaps are not wider than maxGap.
unsigned int fNCell
Number of cells, width of pixel map.