MinervaSegmentationAlg.cxx
Go to the documentation of this file.
1 #include "Geometry/ChannelMapAlgs/MinervaSegmentationAlg.h"
2 #include "fhiclcpp/ParameterSet.h"
3 #include "messagefacility/MessageLogger/MessageLogger.h"
4 #include "Geometry/GeometryCore.h"
5 
6 #include <algorithm>
7 #include <sstream>
8 #include <iostream>
9 #include <cmath>
10 
11 namespace gar {
12  namespace geo {
13  namespace seg {
14 
15  //----------------------------------------------------------------------------
16  /// Default constructor used by derived classes passing the encoding string
17  MinervaSegmentationAlg::MinervaSegmentationAlg(fhicl::ParameterSet const& pset)
18  : SegmentationAlg(pset)
19  {
20  _type = "StripXY";
21  _description = "Cartesian segmentation in the local XY-plane, containing integer number of strips";
22 
23  std::cout << "######### gar::geo::seg::MinervaSegmentationAlg() " << std::endl ;
24 
25  this->reconfigure(pset);
26  }
27 
28  //----------------------------------------------------------------------------
29  /// Default constructor used by derived classes passing an existing decoder
30  MinervaSegmentationAlg::MinervaSegmentationAlg(const BitFieldCoder* decode, fhicl::ParameterSet const& pset)
31  : SegmentationAlg(decode, pset)
32  {
33  _type = "StripXY";
34  _description = "Cartesian segmentation in the local XY-plane, containing integer number of strips";
35 
36  std::cout << "######### gar::geo::seg::MinervaSegmentationAlg() " << std::endl ;
37 
38  this->reconfigure(pset);
39  }
40 
41  //----------------------------------------------------------------------------
42  MinervaSegmentationAlg::~MinervaSegmentationAlg()
43  {
44  }
45 
46  //----------------------------------------------------------------------------
47  void MinervaSegmentationAlg::reconfigure(fhicl::ParameterSet const& pset)
48  {
49  _stripSizeX = pset.get<double>("strip_size_x");
50  _stripSizeY = pset.get<double>("strip_size_y");
51  _encoding = pset.get<std::string>("cellEncoding");
52 
53  _xId = pset.get<std::string>("identifier_x");
54  _yId = pset.get<std::string>("identifier_y");
55  _zId = pset.get<std::string>("identifier_z");
56  _layerId = pset.get<std::string>("identifier_layer");
57  _sliceId = pset.get<std::string>("identifier_slice");
58  _nPlanes = pset.get<unsigned int>("nplanes");
59 
60  _frac = 1./3.;
61 
62  this->PrintParameters();
63 
64  return;
65  }
66 
67  //----------------------------------------------------------------------------
68  void MinervaSegmentationAlg::Initialize(const gar::geo::GeometryCore& )
69  {
70 
71  }
72 
73  //----------------------------------------------------------------------------
74  std::array<double, 3> MinervaSegmentationAlg::GetPosition(const gar::geo::GeometryCore& , const gar::raw::CellID_t& ) const
75  {
76  //Local origin for the Barrel in the middle of the layer
77  //Local origin for the Endcal at the corner of the full stave
78  std::array<double, 3> cellPosition;
79 
80  /* NO OP */
81 
82  return cellPosition;
83  }
84 
85  //----------------------------------------------------------------------------
86  /// determine the cell ID based on the position
87  gar::raw::CellID_t MinervaSegmentationAlg::GetCellID(const gar::geo::GeometryCore& , const unsigned int& det_id, const unsigned int& , const unsigned int& , const unsigned int& layer, const unsigned int& slice, const std::array<double, 3>& localPosition) const
88  {
89  gar::raw::CellID_t cID = 0;
90 
91  _decoder->set(cID, "system", det_id);
92  _decoder->set(cID, "layer", layer);
93  _decoder->set(cID, "slice", slice);
94 
95  double localX = localPosition[0];
96  double localY = localPosition[1];
97  double localZ = localPosition[2];
98 
99  if( localZ < 0 )
100  {
101  //Segmentation in Y
102  //Need to check in which half of the square is the hit --> {y, z} point is below or over the diagonal
103  //(y = z (odd) or y = -z (even))
104 
105  int nCellsX = 1;
106  int nCellsY = int(_layer_dim_Y / (_stripSizeY * 2));
107 
108  int _cellIndexX = int ( localX / ( _layer_dim_X / nCellsX ) );
109  int _cellIndexY = int ( localY / ( _layer_dim_Y / nCellsY ) );
110 
111  //Transform the localX/Y to the local of this cell
112  // double cellOriginX = 0.;
113  double cellOriginY = ( _cellIndexY + 0.5 ) * (_stripSizeY * 2);
114  double cellOriginZ = - 1.;
115 
116  double localy = localY - cellOriginY; //transform it
117  double localz = localZ - cellOriginZ; //transform it
118  bool above = false;
119  bool below = false;
120  if(localy < 0) {
121  //Need to check if the point is below or above y = z
122  if(localy > localz) {
123  above = true;
124  }
125  if(localy < localz) {
126  below = true;
127  }
128  if(above) {
129  _decoder->set(cID, _xId, _cellIndexX);
130  _decoder->set(cID, _yId, _cellIndexY);
131  _decoder->set(cID, _zId, 0);
132  _decoder->set(cID, "triangle", 1);
133  }
134  if(below) {
135  _decoder->set(cID, _xId, _cellIndexX);
136  _decoder->set(cID, _yId, _cellIndexY);
137  _decoder->set(cID, _zId, 0);
138  _decoder->set(cID, "triangle", 0);
139  }
140  }
141  else if( localy > 0 ) {
142  //Need to check if the point is below or above y = -z
143  if(localy > -localz) {
144  above = true;
145  }
146  if(localy < -localz) {
147  below = true;
148  }
149  if(above) {
150  _decoder->set(cID, _xId, _cellIndexX);
151  _decoder->set(cID, _yId, _cellIndexY);
152  _decoder->set(cID, _zId, 0);
153  _decoder->set(cID, "triangle", 3);
154  }
155  if(below) {
156  _decoder->set(cID, _xId, _cellIndexX);
157  _decoder->set(cID, _yId, _cellIndexY);
158  _decoder->set(cID, _zId, 0);
159  _decoder->set(cID, "triangle", 2);
160  }
161  }
162  else {
163  //exception
164  MF_LOG_ERROR("MinervaSegmentationAlg") << "Cannot determine if localX > | < 0";
165  return 0;
166  }
167  } //end if z < 0
168 
169  if( localZ >= 0 )
170  {
171  //Segmentation in X
172  int nCellsX = int(_layer_dim_X / _stripSizeX);
173  int nCellsY = 1;
174 
175  int _cellIndexX = int ( localX / ( _layer_dim_X / nCellsX ) );
176  int _cellIndexY = int ( localY / ( _layer_dim_Y / nCellsY ) );
177 
178  //Transform the localX/Y to the local of this cell
179  double cellOriginX = ( _cellIndexX + 0.5 ) * (_stripSizeX * 2);
180  // double cellOriginY = 0.;
181  double cellOriginZ = 1.;
182 
183  double localx = localX - cellOriginX; //transform it
184  double localz = localZ - cellOriginZ; //transform it
185  bool above = false;
186  bool below = false;
187  if(localx < 0) {
188  //Need to check if the point is below or above y = z
189  if(localx > localz) {
190  above = true;
191  }
192  if(localx < localz) {
193  below = true;
194  }
195  if(above) {
196  _decoder->set(cID, _xId, _cellIndexX);
197  _decoder->set(cID, _yId, _cellIndexY);
198  _decoder->set(cID, _zId, 1);
199  _decoder->set(cID, "triangle", 1);
200  }
201  if(below) {
202  _decoder->set(cID, _xId, _cellIndexX);
203  _decoder->set(cID, _yId, _cellIndexY);
204  _decoder->set(cID, _zId, 1);
205  _decoder->set(cID, "triangle", 0);
206  }
207  }
208  else if( localx > 0 ) {
209  //Need to check if the point is below or above y = -z
210  if(localx > -localz) {
211  above = true;
212  }
213  if(localx < -localz) {
214  below = true;
215  }
216  if(above) {
217  _decoder->set(cID, _xId, _cellIndexX);
218  _decoder->set(cID, _yId, _cellIndexY);
219  _decoder->set(cID, _zId, 1);
220  _decoder->set(cID, "triangle", 3);
221  }
222  if(below) {
223  _decoder->set(cID, _xId, _cellIndexX);
224  _decoder->set(cID, _yId, _cellIndexY);
225  _decoder->set(cID, _zId, 1);
226  _decoder->set(cID, "triangle", 2);
227  }
228  }
229  else {
230  //exception
231  MF_LOG_ERROR("MinervaSegmentationAlg") << "Cannot determine if localX > | < 0";
232  return 0;
233  }
234  } //end if z >= 0
235 
236  return cID;
237  }
238 
239  //------------------------------------------------------------------------------
240  gar::raw::CellID_t MinervaSegmentationAlg::GetComplementaryCellID(gar::raw::CellID_t cellID, unsigned int comp) const
241  {
242  gar::raw::CellID_t cID = 0;
243 
244  if(comp == 0) {
245  //Need to get the lower cellX/Y
246  if(_decoder->get(cellID, "cellX") == 0){
247  _decoder->set(cID, "system", _decoder->get(cellID, "system"));
248  _decoder->set(cID, "layer", _decoder->get(cellID, "layer"));
249  _decoder->set(cID, "slice", _decoder->get(cellID, "slice"));
250  _decoder->set(cID, "cellX", _decoder->get(cellID, "cellX"));
251  _decoder->set(cID, "cellY", _decoder->get(cellID, "cellY")+1);
252  _decoder->set(cID, "cellZ", _decoder->get(cellID, "cellZ"));
253  _decoder->set(cID, "triangle", comp);
254  }
255  if(_decoder->get(cellID, "cellY") == 0){
256  _decoder->set(cID, "system", _decoder->get(cellID, "system"));
257  _decoder->set(cID, "layer", _decoder->get(cellID, "layer"));
258  _decoder->set(cID, "slice", _decoder->get(cellID, "slice"));
259  _decoder->set(cID, "cellX", _decoder->get(cellID, "cellX")+1);
260  _decoder->set(cID, "cellY", _decoder->get(cellID, "cellY"));
261  _decoder->set(cID, "cellZ", _decoder->get(cellID, "cellZ"));
262  _decoder->set(cID, "triangle", comp);
263  }
264  }
265 
266  if(comp == 1 || comp == 2) {
267  //Keep the same cellX/cellY
268  _decoder->set(cID, "system", _decoder->get(cellID, "system"));
269  _decoder->set(cID, "layer", _decoder->get(cellID, "layer"));
270  _decoder->set(cID, "slice", _decoder->get(cellID, "slice"));
271  _decoder->set(cID, "cellX", _decoder->get(cellID, "cellX"));
272  _decoder->set(cID, "cellY", _decoder->get(cellID, "cellY"));
273  _decoder->set(cID, "cellZ", _decoder->get(cellID, "cellZ"));
274  _decoder->set(cID, "triangle", comp);
275  }
276 
277  if(comp == 3) {
278  //Need to get the upper cellX/Y
279  if(_decoder->get(cellID, "cellX") == 0){
280  _decoder->set(cID, "system", _decoder->get(cellID, "system"));
281  _decoder->set(cID, "layer", _decoder->get(cellID, "layer"));
282  _decoder->set(cID, "slice", _decoder->get(cellID, "slice"));
283  _decoder->set(cID, "cellX", _decoder->get(cellID, "cellX"));
284  _decoder->set(cID, "cellY", _decoder->get(cellID, "cellY")-1);
285  _decoder->set(cID, "cellZ", _decoder->get(cellID, "cellZ"));
286  _decoder->set(cID, "triangle", comp);
287  }
288  if(_decoder->get(cellID, "cellY") == 0){
289  _decoder->set(cID, "system", _decoder->get(cellID, "system"));
290  _decoder->set(cID, "layer", _decoder->get(cellID, "layer"));
291  _decoder->set(cID, "slice", _decoder->get(cellID, "slice"));
292  _decoder->set(cID, "cellX", _decoder->get(cellID, "cellX")-1);
293  _decoder->set(cID, "cellY", _decoder->get(cellID, "cellY"));
294  _decoder->set(cID, "cellZ", _decoder->get(cellID, "cellZ"));
295  _decoder->set(cID, "triangle", comp);
296  }
297  }
298 
299  return cID;
300  }
301 
302  //------------------------------------------------------------------------------
303  void MinervaSegmentationAlg::AddHitsMinerva(std::map< gar::raw::CellID_t, std::vector<gar::sdp::CaloDeposit> > &m_Deposits, std::vector<gar::sdp::CaloDeposit> &fDeposits) const
304  {
305  //Loop over the hits in the map and add them together
306  for(auto &it : m_Deposits)
307  {
308  gar::raw::CellID_t cellID = it.first;
309 
310  //Check this cellID
311  //if it is a triangle = 0 or 3
312  //if it is triangle = 1 or 2
313  //need to look for the complementary cellID, add them to this hit and remove from the deposits (to avoid double hit creation)
314 
315  int triangleNb = _decoder->get(cellID, "triangle");
316 
317  //Case bottom or upper triangle, look for complementary cellID
318  if(triangleNb == 0 || triangleNb == 3)
319  {
320  if(triangleNb == 0) {
321  std::vector<gar::sdp::CaloDeposit> vechit = it.second;
322  std::sort(vechit.begin(), vechit.end()); //sort per time
323 
324  float esum = 0.;
325  float stepLength = 0.;
326  float time = vechit.at(0).Time();
327  int trackID = vechit.at(0).TrackID();
328  double pos[3] = { vechit.at(0).X(), vechit.at(0).Y(), vechit.at(0).Z() };
329 
330  for(auto const &hit : vechit) {
331  esum += hit.Energy();
332  stepLength += hit.StepLength();
333  }
334 
335  //need to check if cellX or cellY are - 1
336  gar::raw::CellID_t complementary_cellID = this->GetComplementaryCellID(cellID, 3);
337  auto find = m_Deposits.find( complementary_cellID );
338  if(find != m_Deposits.end()) {
339  std::vector<gar::sdp::CaloDeposit> vechit_comp = find->second;
340  std::sort(vechit_comp.begin(), vechit_comp.end()); //sort per time
341 
342  for(auto const &hit : vechit_comp) {
343  esum += hit.Energy();
344  stepLength += hit.StepLength();
345  }
346 
347  //remove the element from the map now to avoid double counting
348  m_Deposits.erase(find->first);
349  }
350 
351  fDeposits.emplace_back( trackID, time, esum, pos, cellID, stepLength );
352  }
353 
354  if(triangleNb == 3) {
355 
356  std::vector<gar::sdp::CaloDeposit> vechit = it.second;
357  std::sort(vechit.begin(), vechit.end()); //sort per time
358 
359  float esum = 0.;
360  float stepLength = 0.;
361  float time = vechit.at(0).Time();
362  int trackID = vechit.at(0).TrackID();
363  double pos[3] = { vechit.at(0).X(), vechit.at(0).Y(), vechit.at(0).Z() };
364 
365  for(auto const &hit : vechit) {
366  esum += hit.Energy();
367  stepLength += hit.StepLength();
368  }
369 
370  //need to check if cellX or cellY are + 1
371  gar::raw::CellID_t complementary_cellID = this->GetComplementaryCellID(cellID, 0);
372  auto find = m_Deposits.find( complementary_cellID );
373  if(find != m_Deposits.end()) {
374  std::vector<gar::sdp::CaloDeposit> vechit_comp = find->second;
375  std::sort(vechit_comp.begin(), vechit_comp.end()); //sort per time
376 
377  for(auto const &hit : vechit_comp) {
378  esum += hit.Energy();
379  stepLength += hit.StepLength();
380  }
381 
382  //remove the element from the map now to avoid double counting
383  m_Deposits.erase(find->first);
384  }
385 
386  fDeposits.emplace_back( trackID, time, esum, pos, cellID, stepLength );
387  }
388  }
389 
390  if(triangleNb == 1 || triangleNb == 2)
391  {
392  if(triangleNb == 1) {
393  std::vector<gar::sdp::CaloDeposit> vechit = it.second;
394  std::sort(vechit.begin(), vechit.end()); //sort per time
395 
396  float esum = 0.;
397  float stepLength = 0.;
398  float time = vechit.at(0).Time();
399  int trackID = vechit.at(0).TrackID();
400  double pos[3] = { vechit.at(0).X(), vechit.at(0).Y(), vechit.at(0).Z() };
401 
402  for(auto const &hit : vechit) {
403  esum += hit.Energy();
404  stepLength += hit.StepLength();
405  }
406 
407  //need to check if cellX or cellY are the same
408  gar::raw::CellID_t complementary_cellID = this->GetComplementaryCellID(cellID, 2);
409  auto find = m_Deposits.find( complementary_cellID );
410  if(find != m_Deposits.end()) {
411  std::vector<gar::sdp::CaloDeposit> vechit_comp = find->second;
412  std::sort(vechit_comp.begin(), vechit_comp.end()); //sort per time
413 
414  for(auto const &hit : vechit_comp) {
415  esum += hit.Energy();
416  stepLength += hit.StepLength();
417  }
418 
419  //remove the element from the map now to avoid double counting
420  m_Deposits.erase(find->first);
421  }
422  fDeposits.emplace_back( trackID, time, esum, pos, cellID, stepLength );
423  }
424 
425  if(triangleNb == 2) {
426  std::vector<gar::sdp::CaloDeposit> vechit = it.second;
427  std::sort(vechit.begin(), vechit.end()); //sort per time
428 
429  float esum = 0.;
430  float stepLength = 0.;
431  float time = vechit.at(0).Time();
432  int trackID = vechit.at(0).TrackID();
433  double pos[3] = { vechit.at(0).X(), vechit.at(0).Y(), vechit.at(0).Z() };
434 
435  for(auto const &hit : vechit) {
436  esum += hit.Energy();
437  stepLength += hit.StepLength();
438  }
439 
440  //need to check if cellX or cellY are the same
441  gar::raw::CellID_t complementary_cellID = this->GetComplementaryCellID(cellID, 1);
442  auto find = m_Deposits.find( complementary_cellID );
443  if(find != m_Deposits.end()) {
444  std::vector<gar::sdp::CaloDeposit> vechit_comp = find->second;
445  std::sort(vechit_comp.begin(), vechit_comp.end()); //sort per time
446 
447  for(auto const &hit : vechit_comp) {
448  esum += hit.Energy();
449  stepLength += hit.StepLength();
450  }
451 
452  //remove the element from the map now to avoid double counting
453  m_Deposits.erase(find->first);
454  }
455  fDeposits.emplace_back( trackID, time, esum, pos, cellID, stepLength );
456  }
457  }
458 
459  //remove the element from the map now
460  m_Deposits.erase(it.first);
461  }
462  }
463 
464  //----------------------------------------------------------------------------
465  void MinervaSegmentationAlg::PrintParameters() const
466  {
467  std::cout << "cell encoding: " << _encoding << std::endl;
468  std::cout << "identifier_x: " << _xId << std::endl;
469  std::cout << "identifier_y: " << _yId << std::endl;
470  std::cout << "identifier_z: " << _zId << std::endl;
471  std::cout << "strip_size_x: " << _stripSizeX << " cm" << std::endl;
472  std::cout << "strip_size_y: " << _stripSizeY << " cm" << std::endl;
473 
474  std::cout << "identifier_layer: " << _layerId << std::endl;
475  std::cout << "identifier_slice: " << _sliceId << std::endl;
476  std::cout << "nPlanes: " << _nPlanes << std::endl;
477  }
478 
479  }//seg
480  } // geo
481 } //gar
gar::geo::seg::SegmentationAlg::_decoder
const BitFieldCoder * _decoder
Definition: SegmentationAlg.h:126
gar::geo::seg::MinervaSegmentationAlg::MinervaSegmentationAlg
MinervaSegmentationAlg(fhicl::ParameterSet const &pset)
Default constructor used by derived classes passing the encoding string.
Definition: MinervaSegmentationAlg.cxx:17
gar::geo::seg::MinervaSegmentationAlg::_sliceId
std::string _sliceId
the field name used for slice
Definition: MinervaSegmentationAlg.h:95
gar::geo::seg::MinervaSegmentationAlg::_encoding
std::string _encoding
the encoding string
Definition: MinervaSegmentationAlg.h:97
gar::geo::BitFieldCoder::set
void set(long64 &bitfield, size_t index, ulong64 value) const
Definition: BitFieldCoder.h:189
string
std::string string
Definition: nybbler.cc:12
gar::geo::seg::MinervaSegmentationAlg::_yId
std::string _yId
the field name used for Y
Definition: MinervaSegmentationAlg.h:89
gar::geo::seg::MinervaSegmentationAlg::_frac
double _frac
fraction of tiles to remove at the edge
Definition: MinervaSegmentationAlg.h:103
MF_LOG_ERROR
#define MF_LOG_ERROR(category)
Definition: MessageLogger.h:229
gar::geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:436
gar::geo::seg::MinervaSegmentationAlg::AddHitsMinerva
void AddHitsMinerva(std::map< gar::raw::CellID_t, std::vector< gar::sdp::CaloDeposit > > &m_Deposits, std::vector< gar::sdp::CaloDeposit > &fDeposits) const
Definition: MinervaSegmentationAlg.cxx:303
gar::geo::seg::MinervaSegmentationAlg::GetCellID
gar::raw::CellID_t GetCellID(const gar::geo::GeometryCore &geo, const unsigned int &det_id, const unsigned int &stave, const unsigned int &module, const unsigned int &layer, const unsigned int &slice, const std::array< double, 3 > &localPosition) const override
determine the cell ID based on the position
Definition: MinervaSegmentationAlg.cxx:87
gar::geo::seg::MinervaSegmentationAlg::GetComplementaryCellID
gar::raw::CellID_t GetComplementaryCellID(gar::raw::CellID_t cellID, unsigned int comp) const
Definition: MinervaSegmentationAlg.cxx:240
gar::geo::BitFieldCoder
Helper class for decoding and encoding a bit field of 64bits for convenient declaration.
Definition: BitFieldCoder.h:130
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:271
gar::raw::CellID_t
long long int CellID_t
Definition: CaloRawDigit.h:24
gar::geo::seg::MinervaSegmentationAlg::reconfigure
void reconfigure(fhicl::ParameterSet const &pset) override
Definition: MinervaSegmentationAlg.cxx:47
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
gar::geo::seg::MinervaSegmentationAlg::_zId
std::string _zId
the field name used for Z
Definition: MinervaSegmentationAlg.h:91
gar
General GArSoft Utilities.
Definition: anatest_module.cc:57
gar::geo::BitFieldCoder::get
long64 get(long64 bitfield, size_t index) const
Definition: BitFieldCoder.h:181
gar::geo::seg::MinervaSegmentationAlg::GetPosition
std::array< double, 3 > GetPosition(const gar::geo::GeometryCore &geo, const gar::raw::CellID_t &cID) const override
Definition: MinervaSegmentationAlg.cxx:74
gar::geo::seg::MinervaSegmentationAlg::_xId
std::string _xId
the field name used for X
Definition: MinervaSegmentationAlg.h:87
CLHEP::decode
void decode(std::any const &a, Hep2Vector &result)
Definition: CLHEP_ps.h:12
geo
LArSoft geometry interface.
Definition: ChannelGeo.h:16
gar::geo::seg::MinervaSegmentationAlg::_layerId
std::string _layerId
the field name used for layer
Definition: MinervaSegmentationAlg.h:93
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
gar::geo::seg::MinervaSegmentationAlg::_nPlanes
unsigned int _nPlanes
number of planes
Definition: MinervaSegmentationAlg.h:105
gar::geo::seg::MinervaSegmentationAlg::Initialize
void Initialize(const gar::geo::GeometryCore &geo) override
Definition: MinervaSegmentationAlg.cxx:68
Generated by   doxygen 1.8.11