Public Member Functions | Private Types | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
gar::rec::alg::StripSplitterAlg Class Reference

#include <StripSplitterAlg.h>

Public Member Functions

 StripSplitterAlg (fhicl::ParameterSet const &pset)
 
virtual ~StripSplitterAlg ()
 
void reconfigure (fhicl::ParameterSet const &pset)
 
void ClearLists ()
 
void PrepareAlgo (const std::vector< art::Ptr< gar::rec::CaloHit > > &hitVector)
 
void DoStripSplitting ()
 
bool GetSaveStripEndsFlag () const
 
std::vector< const gar::rec::CaloHit * > getStripEndsHits () const
 
std::vector< const gar::rec::CaloHit * > getIntersectionHits () const
 
std::vector< const gar::rec::CaloHit * > getSplitHits () const
 
std::vector< const gar::rec::CaloHit * > getUnSplitHits () const
 

Private Types

enum  StripOrientation { TRANSVERSE = 0, LONGITUDINAL }
 

Private Member Functions

void getVirtualHits (const gar::rec::CaloHit *hit, int orientation, bool isBarrel, std::vector< const gar::rec::CaloHit * > &virtualhits)
 
TVector3 stripIntersect (const gar::rec::CaloHit *hit0, const TVector3 &dir0, const gar::rec::CaloHit *hit1, const TVector3 &dir1)
 

Static Private Member Functions

static bool SortByLayer (const gar::rec::CaloHit *rha, const gar::rec::CaloHit *rhb)
 

Private Attributes

gar::geo::GeometryCore const * fGeo
 geometry information More...
 
gar::geo::BitFieldCoder const * fFieldDecoder
 
std::string fSSAAlgName
 
std::string fDet
 
bool fSaveStripEnds
 
int fInnerSymmetry
 
std::string fEncoding
 
double fStripWidth
 
double fStripLength
 
int fnVirtual
 
std::vector< const gar::rec::CaloHit * > m_CaloHitVecOdd
 
std::vector< const gar::rec::CaloHit * > m_CaloHitVecEven
 
std::vector< const gar::rec::CaloHit * > fStripEndsHits
 
std::vector< const gar::rec::CaloHit * > fIntersectionHits
 
std::vector< const gar::rec::CaloHit * > unSplitStripHits
 
std::vector< const gar::rec::CaloHit * > splitStripHits
 

Detailed Description

Definition at line 27 of file StripSplitterAlg.h.

Member Enumeration Documentation

Enumerator
TRANSVERSE 
LONGITUDINAL 

Definition at line 81 of file StripSplitterAlg.h.

Constructor & Destructor Documentation

gar::rec::alg::StripSplitterAlg::StripSplitterAlg ( fhicl::ParameterSet const &  pset)

Definition at line 26 of file StripSplitterAlg.cxx.

27  {
28  ClearLists();
29 
30  fGeo = gar::providerFrom<geo::GeometryGAr>();
31  fInnerSymmetry = 0;
32  fStripWidth = -1; //cm
33  fStripLength = -1;//cm
34  fnVirtual = 1;
35 
36  this->reconfigure(pset);
37 
38  return;
39  }
gar::geo::GeometryCore const * fGeo
geometry information
void reconfigure(fhicl::ParameterSet const &pset)
gar::rec::alg::StripSplitterAlg::~StripSplitterAlg ( )
virtual

Definition at line 42 of file StripSplitterAlg.cxx.

43  {
44  if ( fFieldDecoder ) delete fFieldDecoder;
45  return;
46  }
gar::geo::BitFieldCoder const * fFieldDecoder

Member Function Documentation

void gar::rec::alg::StripSplitterAlg::ClearLists ( )

Definition at line 74 of file StripSplitterAlg.cxx.

75  {
76  m_CaloHitVecOdd.clear();
77  m_CaloHitVecEven.clear();
78 
79  fStripEndsHits.clear();
80  fIntersectionHits.clear();
81 
82  unSplitStripHits.clear();
83  splitStripHits.clear();
84 
85  return;
86  }
std::vector< const gar::rec::CaloHit * > m_CaloHitVecEven
std::vector< const gar::rec::CaloHit * > m_CaloHitVecOdd
std::vector< const gar::rec::CaloHit * > fIntersectionHits
std::vector< const gar::rec::CaloHit * > splitStripHits
std::vector< const gar::rec::CaloHit * > fStripEndsHits
std::vector< const gar::rec::CaloHit * > unSplitStripHits
void gar::rec::alg::StripSplitterAlg::DoStripSplitting ( )

Definition at line 126 of file StripSplitterAlg.cxx.

127  {
128  MF_LOG_DEBUG("StripSplitterAlg")
129  << "StripSplitterAlg::DoStripSplitting()";
130 
131  //Collection to split
132  std::vector <const gar::rec::CaloHit*> *toSplit;
133  int orientation;
134  //Orientation
135  //even layers are segmented in Y (local) -> Transverse
136  //odd layers are segmented in X (local) -> Longitudinal
137 
138  // loop over strip collections in even and odd layers (assumed to have perpendicular orientations)
139  for (int icol = 0; icol < 2; icol++)
140  {
141  MF_LOG_DEBUG("StripSplitterAlg")
142  << "Start Loop " << icol;
143 
144  switch (icol)
145  {
146  case 0:
147  orientation = TRANSVERSE;//even layers
148  toSplit = &m_CaloHitVecEven;
149  break;
150  case 1:
151  orientation = LONGITUDINAL;//odd layers
152  toSplit = &m_CaloHitVecOdd;
153  break;
154  default:
155  MF_LOG_ERROR("StripSplitterAlg::DoStripSplitting")
156  << "Crazy stuff!";
157  throw cet::exception("StripSplitterAlg::DoStripSplitting");
158  }
159 
160  // loop over hits of this type (long/trans) collection
161  for (uint i = 0; i < toSplit->size(); i++)
162  {
163  const gar::rec::CaloHit* hit = toSplit->at(i);
164  //1 == Barrel, 2 == Endcap, 4 == MuID
165  unsigned int det_id = fFieldDecoder->get(hit->CellID(), "system");
166 
167  if ( det_id != 1 && det_id != 2 && det_id != 4 )
168  {
169  MF_LOG_ERROR("StripSplitterAlg::DoStripSplitting")
170  << " Check det it " << det_id
171  << " Problem with Hit " << i
172  << " pointing at " << hit
173  << " at position ( " << hit->Position()[0] << " cm , " << hit->Position()[1] << " cm , " << hit->Position()[2] << " cm )"
174  << " orientation " << (orientation == LONGITUDINAL ? "LONGITUDINAL" : "TRANSVERSE");
175  throw cet::exception("StripSplitterAlg::DoStripSplitting");
176  }
177 
178  MF_LOG_DEBUG("StripSplitterAlg::DoStripSplitting") << "recohit " << hit
179  << " with cellID " << hit->CellID()
180  << " has energy " << hit->Energy() * CLHEP::MeV / CLHEP::GeV
181  << " pos (" << hit->Position()[0] << ", " << hit->Position()[1] << ", " << hit->Position()[2] << ")";
182 
183  // split the hits
184  std::vector <const gar::rec::CaloHit*> virtualhits;
185  bool isBarrel = (det_id == 0 || det_id == 4) ? true : false;
186  /* DEBLEO: Should it be bool isBarrel = (det_id == 1); */
187  getVirtualHits(hit, orientation, isBarrel, virtualhits);
188 
189  // add (new) hits to collections
190  if (virtualhits.size() == 0)
191  {
192 
193  MF_LOG_DEBUG("StripSplitterAlg")
194  << " Adding unsplit hit1 " << i
195  << " pointing at " << hit
196  << " orientation " << (orientation == LONGITUDINAL ? "LONGITUDINAL" : "TRANSVERSE");
197 
198  // not split, add original hit
199  unSplitStripHits.emplace_back(hit);
200  } else {
201  // split was split, add the virtual hits
202  for (uint hh = 0; hh < virtualhits.size(); hh++)
203  {
204 
205  MF_LOG_DEBUG("StripSplitterAlg")
206  << "adding virtual hit " << hh;
207 
208  splitStripHits.emplace_back(virtualhits.at(hh));
209  }
210  }
211  } // loop over hits
212  } // long/trans loop
213 
214  return;
215  }
raw::CellID_t CellID() const
Definition: CaloHit.h:72
std::vector< const gar::rec::CaloHit * > m_CaloHitVecEven
std::vector< const gar::rec::CaloHit * > m_CaloHitVecOdd
#define MF_LOG_ERROR(category)
static constexpr double MeV
Definition: Units.h:129
gar::geo::BitFieldCoder const * fFieldDecoder
static constexpr double GeV
Definition: Units.h:28
float Energy() const
Definition: CaloHit.h:69
std::vector< const gar::rec::CaloHit * > splitStripHits
Detector simulation of raw signals on wires.
long64 get(long64 bitfield, size_t index) const
const float * Position() const
Definition: CaloHit.h:70
#define MF_LOG_DEBUG(id)
std::vector< const gar::rec::CaloHit * > unSplitStripHits
unsigned uint
Definition: qglobal.h:351
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void getVirtualHits(const gar::rec::CaloHit *hit, int orientation, bool isBarrel, std::vector< const gar::rec::CaloHit * > &virtualhits)
std::vector<const gar::rec::CaloHit*> gar::rec::alg::StripSplitterAlg::getIntersectionHits ( ) const
inline

Definition at line 46 of file StripSplitterAlg.h.

46 { return fIntersectionHits; }
std::vector< const gar::rec::CaloHit * > fIntersectionHits
bool gar::rec::alg::StripSplitterAlg::GetSaveStripEndsFlag ( ) const
inline

Definition at line 42 of file StripSplitterAlg.h.

std::vector<const gar::rec::CaloHit*> gar::rec::alg::StripSplitterAlg::getSplitHits ( ) const
inline

Definition at line 48 of file StripSplitterAlg.h.

48 { return splitStripHits; }
std::vector< const gar::rec::CaloHit * > splitStripHits
std::vector<const gar::rec::CaloHit*> gar::rec::alg::StripSplitterAlg::getStripEndsHits ( ) const
inline

Definition at line 44 of file StripSplitterAlg.h.

44 { return fStripEndsHits; }
std::vector< const gar::rec::CaloHit * > fStripEndsHits
std::vector<const gar::rec::CaloHit*> gar::rec::alg::StripSplitterAlg::getUnSplitHits ( ) const
inline

Definition at line 50 of file StripSplitterAlg.h.

50 { return unSplitStripHits; }
std::vector< const gar::rec::CaloHit * > unSplitStripHits
void gar::rec::alg::StripSplitterAlg::getVirtualHits ( const gar::rec::CaloHit hit,
int  orientation,
bool  isBarrel,
std::vector< const gar::rec::CaloHit * > &  virtualhits 
)
private

Definition at line 218 of file StripSplitterAlg.cxx.

221  {
222  MF_LOG_DEBUG("StripSplitterAlg")
223  << "StripSplitterAlg::getVirtualHits()";
224 
225  // this splits the strip into zero or more hits along its length
226  // by looking at nearby hits with different orientation (trans/long or square)
227 
228  int detid = fFieldDecoder->get(hit->CellID(), "system");
229  int layer = hit->Layer();
230  int module = fFieldDecoder->get(hit->CellID(), "module");
231  int stave = fFieldDecoder->get(hit->CellID(), "stave");
232 
233  const std::array<double, 3> pt = { hit->Position()[0], hit->Position()[1], hit->Position()[2] };
235  fStripLength = fGeo->getStripLength(pt, hit->CellID());
236  fnVirtual = 0;
237 
238  if(fStripWidth != 0)
240 
241  MF_LOG_DEBUG("StripSplitterAlg")
242  << " StripSplitterAlg::getVirtualHits()"
243  << " Strip splitted in " << fnVirtual << " virtual cells";
244 
245  // get the ends of this strip
246  std::pair < TVector3, TVector3 > stripEnds = fGeo->GetStripEnds(pt, hit->CellID());
247  TVector3 stripDir = stripEnds.first - stripEnds.second;
248 
249  if(fSaveStripEnds)
250  {
251  float ppp[3];
252  ppp[0] = stripEnds.first.X();
253  ppp[1] = stripEnds.first.Y();
254  ppp[2] = stripEnds.first.Z();
255 
256  const gar::rec::CaloHit* StripEndHit1 = new gar::rec::CaloHit(0.01, 0., ppp, 0., 0);
257  fStripEndsHits.emplace_back(StripEndHit1);
258 
259  ppp[0] = stripEnds.second.X();
260  ppp[1] = stripEnds.second.Y();
261  ppp[2] = stripEnds.second.Z();
262 
263  const gar::rec::CaloHit* StripEndHit2 = new gar::rec::CaloHit(0.01, 0., ppp, 0., 0);
264  fStripEndsHits.emplace_back(StripEndHit2);
265  }
266 
267  // decide which collections to use to split the strip
268  int splitterOrientation;
269  std::vector <const gar::rec::CaloHit*> *splitterVec;
270 
271  if ( layer%2 == 0 )
272  {
273  splitterVec = &m_CaloHitVecOdd;
274  splitterOrientation = LONGITUDINAL;
275  }
276  else
277  {
278  splitterVec = &m_CaloHitVecEven;
279  splitterOrientation = TRANSVERSE;
280  }
281 
282  std::map <int, float> virtEnergy;
283  int nSplitters(0);
284 
285  // loop over splitter cols, find nearby hits
286  // strips, cells
287  for (int jj = 0; jj < 2; jj++)
288  {
289  std::vector <const gar::rec::CaloHit*> *splitter = splitterVec;
290 
291  for (uint i = 0; i < splitter->size(); i++)
292  {
293  const gar::rec::CaloHit *hit2 = splitter->at(i);
294 
295  int detid2 = fFieldDecoder->get(hit2->CellID(), "system");
296  int layer2 = hit2->Layer();
297  int module2 = fFieldDecoder->get(hit2->CellID(), "module");
298  int stave2 = fFieldDecoder->get(hit2->CellID(), "stave");
299 
300  int ddetid = std::abs(detid2 - detid);
301  int dlayer = std::abs(layer2 - layer);
302  int dstave = std::abs(stave2 - stave);
303  int dmodule = std::abs(module2 - module);
304 
305  TVector3 point2(hit2->Position()[0], hit2->Position()[1], hit2->Position()[2]);
306 
307  //if they are in different part of the det, ignore
308  if(ddetid != 0) continue;
309 
310  //TODO CHECK!
311  // are the two hits close enough to look at further?
312  // if hits in same module and same stave, require that only one layer difference
313  if (dmodule == 0 && dstave == 0 && dlayer > 1) continue;
314 
315  if (isBarrel)
316  {
317  dstave = std::min( dstave, fInnerSymmetry - dstave);
318  if ( dstave == 0 && dmodule > 1 ) continue; // allow same stave and +- 1 module
319  if ( dmodule == 0 && dstave > 1 ) continue; // or same module +- 1 stave
320  if ( dstave == 0 && dlayer > 1) continue; // if in same stave, require dlayer==1
321  }
322  else
323  {
324  //module == 0 or 6 for endcap
325  //stave is from 1 to 4
326  //TODO Check this!
327  //endcap
328  dstave = std::min( dstave, 4 - dstave);
329  if (dmodule != 0) continue; // different endcap
330  if (dstave > 1) continue; // more than 1 stave (=quarter endcap) apart
331  if (dlayer > 1) continue; // more than 1 layer apart
332  }
333 
334  MF_LOG_DEBUG("StripSplitterAlg::getVirtualHits()")
335  << " Getting hit2 " << i
336  << " pointing at " << hit2
337  << " orientation " << (splitterOrientation == LONGITUDINAL ? "LONGITUDINAL" : "TRANSVERSE")
338  << " dtave " << dstave
339  << " dmodule " << dmodule
340  << " dlayer " << dlayer;
341 
342  // simple distance check for remaining hit pairs
343  float dist = std::sqrt( std::pow(hit2->Position()[0] - hit->Position()[0], 2) + std::pow(hit2->Position()[1] - hit->Position()[1], 2) + std::pow(hit2->Position()[2] - hit->Position()[2], 2) );
344 
345  if (dist > 2*fStripLength) {
346  MF_LOG_DEBUG("StripSplitterAlg::getVirtualHits()")
347  << " Distance between hit1 and hit2 " << dist
348  << " > 2*fStripLength " << 2*fStripLength;
349  continue;
350  }
351 
352  // for remaining hits, check if they overlap
353  TVector3 stripDir2(0, 0, 0);
354  if (jj == 0) {
355  //strip
356  std::array<double, 3> pt2 = { hit2->Position()[0], hit2->Position()[1], hit2->Position()[2] };
357  std::pair < TVector3, TVector3 > stripEnds2 = fGeo->GetStripEnds(pt2, hit2->CellID());
358  stripDir2 = stripEnds2.first - stripEnds2.second;
359  }
360 
361  // check if strips intersect
362  TVector3 intercept = stripIntersect(hit, stripDir, hit2, stripDir2);
363  // intercept found, calculate in which virtual cell
364  if (intercept.Mag() > 0)
365  {
366  nSplitters++;
367  float frac(-1);
368  for (int ii = 0; ii < 3; ii++) {
369  float dx = stripEnds.second[ii] - stripEnds.first[ii];
370  if (std::fabs(dx) > 0.1) {
371  frac = ( intercept[ii] - stripEnds.first[ii] ) / dx;
372  break;
373  }
374  }
375 
376  if (frac >= 0.0 && frac <= 1.0)
377  {
378  int segment = int(frac * fnVirtual);
379  if (segment >= 0 && segment < fnVirtual)
380  {
381  if (virtEnergy.find(segment) != virtEnergy.end())
382  {
383  virtEnergy[segment] += hit2->Energy();
384  } else {
385  virtEnergy[segment] = hit2->Energy();
386  }
387 
388  if(fSaveStripEnds)
389  {
390  float pos[3];
391  pos[0] = intercept.X();
392  pos[1] = intercept.Y();
393  pos[2] = intercept.Z();
394 
395  const gar::rec::CaloHit* interhit = new gar::rec::CaloHit(0.1, 0., pos, 0., 0);
396  fIntersectionHits.emplace_back(interhit);
397  }
398 
399  } else {
400  MF_LOG_DEBUG ("StripSplitterAlg::getVirtualHits()")
401  << "strange segment " << segment
402  << " frac = " << frac
403  << " nvirt = " << fnVirtual;
404  }
405  } else {
406  MF_LOG_DEBUG ("StripSplitterAlg::getVirtualHits()")
407  << "strange frac " << frac;
408  }
409  }
410  }
411  }
412 
413  MF_LOG_DEBUG("StripSplitterAlg::getVirtualHits()")
414  << " Number of splitters " << nSplitters;
415 
416  // now create the virtual cells, and assign energy
417  float totenergy(0);
418  for (std::map <int, float>::iterator it = virtEnergy.begin(); it != virtEnergy.end(); ++it) {
419  totenergy += it->second;
420  }
421 
422  for (std::map <int, float>::iterator it = virtEnergy.begin(); it != virtEnergy.end(); ++it) {
423  // energy of hit
424  float energy = hit->Energy() * it->second / totenergy;
425  // position of hit
426  TVector3 virtualCentre = stripEnds.second - stripEnds.first;
427  virtualCentre *= (it->first + 0.5) / fnVirtual;
428  virtualCentre += stripEnds.first;
429 
430  float pos[3];
431  pos[0] = virtualCentre.X();
432  pos[1] = virtualCentre.Y();
433  pos[2] = virtualCentre.Z();
434 
435  // make the new hit
436  const gar::rec::CaloHit* newhit = new gar::rec::CaloHit(energy, hit->Time(), pos, hit->CellID(), hit->Layer());
437 
438  MF_LOG_DEBUG("StripSplitterAlg::getVirtualHits()")
439  << " Creating new virtual hit pointing at " << newhit
440  << " Energy " << energy
441  << " Position " << pos[0] << " " << pos[1] << " " << pos[2];
442 
443  virtualhits.emplace_back(newhit);
444  }
445  }
gar::geo::GeometryCore const * fGeo
geometry information
intermediate_table::iterator iterator
std::pair< TVector3, TVector3 > GetStripEnds(const std::array< double, 3 > &point, const gar::raw::CellID_t &cID) const
raw::CellID_t CellID() const
Definition: CaloHit.h:72
TVector3 stripIntersect(const gar::rec::CaloHit *hit0, const TVector3 &dir0, const gar::rec::CaloHit *hit1, const TVector3 &dir1)
std::vector< const gar::rec::CaloHit * > m_CaloHitVecEven
std::vector< const gar::rec::CaloHit * > m_CaloHitVecOdd
double getStripWidth(const std::array< double, 3 > &point) const
std::pair< float, float > Time() const
Definition: CaloHit.h:71
constexpr T pow(T x)
Definition: pow.h:72
double getStripLength(const std::array< double, 3 > &point, const gar::raw::CellID_t &cID) const
std::vector< const gar::rec::CaloHit * > fIntersectionHits
unsigned int Layer() const
Definition: CaloHit.h:73
T abs(T value)
gar::geo::BitFieldCoder const * fFieldDecoder
float Energy() const
Definition: CaloHit.h:69
std::vector< const gar::rec::CaloHit * > fStripEndsHits
long64 get(long64 bitfield, size_t index) const
const float * Position() const
Definition: CaloHit.h:70
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
#define MF_LOG_DEBUG(id)
unsigned uint
Definition: qglobal.h:351
void gar::rec::alg::StripSplitterAlg::PrepareAlgo ( const std::vector< art::Ptr< gar::rec::CaloHit > > &  hitVector)

Definition at line 95 of file StripSplitterAlg.cxx.

96  {
97  MF_LOG_DEBUG("StripSplitterAlg")
98  << "StripSplitterAlg::PrepareAlgo()";
99 
100  //Clear the lists
101  ClearLists();
102 
103  //Loop over all hits
104  for (std::vector< art::Ptr<gar::rec::CaloHit> >::const_iterator iter = hitVector.begin(), iterEnd = hitVector.end(); iter != iterEnd; ++iter)
105  {
106  art::Ptr<gar::rec::CaloHit> hitPtr = *iter;
107  const gar::rec::CaloHit *hit = hitPtr.get();
108 
109  //check the layer
110  unsigned int layer = hit->Layer();
111  if(layer%2 == 0)
112  m_CaloHitVecEven.emplace_back( hit );
113  else{
114  m_CaloHitVecOdd.emplace_back( hit );
115  }
116  }
117 
118  //Sort the hits by layer
121 
122  return;
123  }
std::vector< const gar::rec::CaloHit * > m_CaloHitVecEven
std::vector< const gar::rec::CaloHit * > m_CaloHitVecOdd
struct vector vector
static bool SortByLayer(const gar::rec::CaloHit *rha, const gar::rec::CaloHit *rhb)
unsigned int Layer() const
Definition: CaloHit.h:73
Detector simulation of raw signals on wires.
#define MF_LOG_DEBUG(id)
T const * get() const
Definition: Ptr.h:149
void gar::rec::alg::StripSplitterAlg::reconfigure ( fhicl::ParameterSet const &  pset)

Definition at line 49 of file StripSplitterAlg.cxx.

50  {
51  fSSAAlgName = pset.get<std::string>("SSAAlgName");
52  fDet = pset.get<std::string>("DetectorSystem"); //ECAL or MuID
53  fSaveStripEnds = pset.get<bool>("SaveStripEnds", false);
54 
55  if(fDet == "ECAL") {
59  }
60  else if(fDet == "MuID") {
64  } else {
65  MF_LOG_ERROR("StripSplitterAlg::reconfigure")
66  << "Detector system parameter incorrectly set. Should be ECAL or MuID -- Exiting!";
67  throw cet::exception("StripSplitterAlg::reconfigure");
68  }
69 
70  return;
71  }
gar::geo::GeometryCore const * fGeo
geometry information
const std::string GetMuIDCellIDEncoding() const
int GetECALInnerSymmetry() const
Definition: GeometryCore.h:957
std::string string
Definition: nybbler.cc:12
#define MF_LOG_ERROR(category)
gar::geo::BitFieldCoder const * fFieldDecoder
Helper class for decoding and encoding a bit field of 64bits for convenient declaration.
const std::string GetECALCellIDEncoding() const
int GetMuIDInnerSymmetry() const
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool gar::rec::alg::StripSplitterAlg::SortByLayer ( const gar::rec::CaloHit rha,
const gar::rec::CaloHit rhb 
)
staticprivate

Definition at line 89 of file StripSplitterAlg.cxx.

90  {
91  return rha->Layer() < rhb->Layer() ;
92  }
unsigned int Layer() const
Definition: CaloHit.h:73
TVector3 gar::rec::alg::StripSplitterAlg::stripIntersect ( const gar::rec::CaloHit hit0,
const TVector3 &  dir0,
const gar::rec::CaloHit hit1,
const TVector3 &  dir1 
)
private

Definition at line 448 of file StripSplitterAlg.cxx.

449  {
450  // find intercept of hit1 with hit0
451  // hit0 must be a strip
452  // hit1 can be an orthogonal strip, or a cell
453  // dir0,1 are direction of strip
454 
455  // centre position of cell/strip
456  TVector3 stripCentre[2];
457  stripCentre[0].SetXYZ( hit0->Position()[0], hit0->Position()[1], hit0->Position()[2] );
458  stripCentre[1].SetXYZ( hit1->Position()[0], hit1->Position()[1], hit1->Position()[2] );
459 
460  const std::array<double, 3> pt = { hit0->Position()[0], hit0->Position()[1], hit0->Position()[2] };
461  fStripLength = fGeo->getStripLength(pt, hit0->CellID());
462 
463  // direction of strip long axis
464  // 0,0,0 for square cell
465  TVector3 stripDir[2];
466  stripDir[0] = dir0;
467  stripDir[1] = dir1;
468 
469  // deal with cell case
470  // define it's direction as perpendicular to strip and vector cell centre to origin
471  bool isStrip[2];
472  for (int i = 0; i < 2; i++)
473  {
474  if ( stripDir[i].Mag() > 1e-10 ) {
475  isStrip[i] = true;
476  } else {
477  isStrip[i] = false;
478  stripDir[i] = stripDir[1-i].Cross(stripCentre[i]);
479  }
480  // ensure dir is normalised
481  stripDir[i] *= 1. / stripDir[i].Mag();
482  }
483 
484  if (not isStrip[0]) {
485  MF_LOG_ERROR ("StripSplitterAlg::stripIntersect")
486  << "first hit should be a strip";
487  throw cet::exception("StripSplitterAlg::stripIntersect");
488  }
489 
490  TVector3 p[2][2]; // ends of strips
491  for (int j = 0; j < 2; j++) {
492  for (int i = 0; i < 2; i++) {
493  float ll = isStrip[j] ? fStripLength : fStripWidth*1.1; // inflate a little for cells
494  p[j][i] = stripCentre[j] - std::pow(-1, i) * 0.5 * ll * stripDir[j];
495  }
496  }
497 
498  TVector3 pNorm[2][2];
499  for (int j = 0; j < 2; j++) {
500  for (int i = 0; i < 2; i++) {
501  float mm = p[j][i].Mag();
502  pNorm[j][i] = p[j][i];
503  pNorm[j][i] *= 1./mm;
504  }
505  }
506 
507  TVector3 inPlane[2]; // difference between points: line inside plane
508  for (int j = 0; j < 2; j++) {
509  inPlane[j] = p[j][0] - p[j][1];
510  }
511 
512  // vector normal to both lines (this is normal to the plane we're interested in)
513  TVector3 normal = inPlane[0].Cross(inPlane[1]);
514  float mag = normal.Mag();
515  normal *= 1./mag;
516 
517  // point on line [0]
518  TVector3 point(p[0][0] + p[0][1]);
519  point *= 0.5;
520 
521  // calculate the projected positions of ends of [1] on
522  // the plane which contains "point" with normal "normal"
523  TVector3 qPrime[2];
524  for (int i = 0; i < 2; i++) {
525  float d = ((point - p[1][i]).Dot(normal)) / (pNorm[1][i].Dot(normal));
526  qPrime[i] = p[1][i] + d * pNorm[1][i];
527  }
528 
529  // find the intersection of lines qPrime and p (they are in same plane)
530  TVector3 a(p[0][1] - p[0][0]);
531  TVector3 b(qPrime[1] - qPrime[0]);
532  TVector3 c(qPrime[0] - p[0][0]);
533 
534  float factor = ( c.Cross(b) ).Dot( a.Cross(b) ) / ( ( a.Cross(b) ).Mag2() );
535 
536  TVector3 x = p[0][0] + a*factor;
537 
538  // check that two lines really intercept
539  bool intersect = true;
540  for (int ii = 0; ii < 2; ii++) {
541  for (int j = 0; j < 3; j++) {
542  float d0 = ii ==0 ? x[j] - p[0][0][j] : x[j] - qPrime[0][j];
543  float d1 = ii ==0 ? x[j] - p[0][1][j] : x[j] - qPrime[1][j];
544  if ( d0 * d1 > 1e-10 ) {
545  intersect = false;
546  }
547  }
548  }
549 
550  MF_LOG_DEBUG ("StripSplitterAlg::stripIntersect")
551  << "Intersection found " << intersect;
552 
553  if (intersect) return x;
554  else return TVector3(0,0,0);
555  }
gar::geo::GeometryCore const * fGeo
geometry information
raw::CellID_t CellID() const
Definition: CaloHit.h:72
constexpr T pow(T x)
Definition: pow.h:72
#define MF_LOG_ERROR(category)
double getStripLength(const std::array< double, 3 > &point, const gar::raw::CellID_t &cID) const
const double e
const double a
p
Definition: test.py:223
const float * Position() const
Definition: CaloHit.h:70
#define MF_LOG_DEBUG(id)
static constexpr double mm
Definition: Units.h:65
static bool * b
Definition: config.cpp:1043
list x
Definition: train.py:276
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Member Data Documentation

std::string gar::rec::alg::StripSplitterAlg::fDet
private

Definition at line 64 of file StripSplitterAlg.h.

std::string gar::rec::alg::StripSplitterAlg::fEncoding
private

Definition at line 68 of file StripSplitterAlg.h.

gar::geo::BitFieldCoder const* gar::rec::alg::StripSplitterAlg::fFieldDecoder
private

Definition at line 61 of file StripSplitterAlg.h.

gar::geo::GeometryCore const* gar::rec::alg::StripSplitterAlg::fGeo
private

geometry information

Definition at line 60 of file StripSplitterAlg.h.

int gar::rec::alg::StripSplitterAlg::fInnerSymmetry
private

Definition at line 67 of file StripSplitterAlg.h.

std::vector<const gar::rec::CaloHit*> gar::rec::alg::StripSplitterAlg::fIntersectionHits
private

Definition at line 77 of file StripSplitterAlg.h.

int gar::rec::alg::StripSplitterAlg::fnVirtual
private

Definition at line 71 of file StripSplitterAlg.h.

bool gar::rec::alg::StripSplitterAlg::fSaveStripEnds
private

Definition at line 65 of file StripSplitterAlg.h.

std::string gar::rec::alg::StripSplitterAlg::fSSAAlgName
private

Definition at line 63 of file StripSplitterAlg.h.

std::vector<const gar::rec::CaloHit*> gar::rec::alg::StripSplitterAlg::fStripEndsHits
private

Definition at line 76 of file StripSplitterAlg.h.

double gar::rec::alg::StripSplitterAlg::fStripLength
private

Definition at line 70 of file StripSplitterAlg.h.

double gar::rec::alg::StripSplitterAlg::fStripWidth
private

Definition at line 69 of file StripSplitterAlg.h.

std::vector<const gar::rec::CaloHit*> gar::rec::alg::StripSplitterAlg::m_CaloHitVecEven
private

Definition at line 74 of file StripSplitterAlg.h.

std::vector<const gar::rec::CaloHit*> gar::rec::alg::StripSplitterAlg::m_CaloHitVecOdd
private

Definition at line 73 of file StripSplitterAlg.h.

std::vector<const gar::rec::CaloHit*> gar::rec::alg::StripSplitterAlg::splitStripHits
private

Definition at line 79 of file StripSplitterAlg.h.

std::vector<const gar::rec::CaloHit*> gar::rec::alg::StripSplitterAlg::unSplitStripHits
private

Definition at line 78 of file StripSplitterAlg.h.


The documentation for this class was generated from the following files: