Public Member Functions | Private Member Functions | Private Attributes | List of all members
opdet::OpSlicer Class Reference
Inheritance diagram for opdet::OpSlicer:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Member Functions

 OpSlicer (const fhicl::ParameterSet &)
 
virtual ~OpSlicer ()
 
void beginJob ()
 
void endJob ()
 
void reconfigure (fhicl::ParameterSet const &pset)
 
void produce (art::Event &)
 
- Public Member Functions inherited from art::EDProducer
 EDProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Produceroperator= (Producer const &)=delete
 
Produceroperator= (Producer &&)=delete
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifieroperator= (Modifier const &)=delete
 
Modifieroperator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Member Functions

double YZDist (art::Ptr< recob::OpHit >, art::Ptr< recob::OpHit >) const
 
double Dist (art::Ptr< recob::OpHit >, art::Ptr< recob::OpHit >) const
 
int YZCentroid (std::vector< art::Ptr< recob::OpHit > >, std::vector< int >) const
 
void GetHitYZ (std::vector< art::Ptr< recob::OpHit > >, std::vector< int >, std::vector< double > &, std::vector< double > &) const
 
void ClusterHits (std::vector< art::Ptr< recob::OpHit >>, std::vector< recob::OpFlash > &, std::vector< std::vector< int > > &, detinfo::DetectorClocksData const &) const
 

Private Attributes

std::string fOpHitModuleLabel
 
double fTScale
 
double fRScale
 
double fR0
 
double fBreakTime
 
int fMinN
 
double fTrigCoinc
 
art::ServiceHandle< geo::Geometrygeo
 

Additional Inherited Members

- Public Types inherited from art::EDProducer
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Static Public Member Functions inherited from art::EDProducer
static void commitEvent (EventPrincipal &ep, Event &e)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollectorconsumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 50 of file OpSlicer_module.cc.

Constructor & Destructor Documentation

opdet::OpSlicer::OpSlicer ( const fhicl::ParameterSet pset)
explicit

Definition at line 107 of file OpSlicer_module.cc.

107  : EDProducer{pset}
108  {
109 
110  reconfigure(pset);
111 
112  produces< std::vector< recob::OpFlash > >();
113  produces< art::Assns< recob::OpFlash, recob::OpHit > >();
114 
115  }
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
void reconfigure(fhicl::ParameterSet const &pset)
opdet::OpSlicer::~OpSlicer ( )
virtual

Definition at line 136 of file OpSlicer_module.cc.

137  {
138  }

Member Function Documentation

void opdet::OpSlicer::beginJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 140 of file OpSlicer_module.cc.

141  {
142  }
void opdet::OpSlicer::ClusterHits ( std::vector< art::Ptr< recob::OpHit >>  ,
std::vector< recob::OpFlash > &  ,
std::vector< std::vector< int > > &  ,
detinfo::DetectorClocksData const &   
) const
private

Definition at line 262 of file OpSlicer_module.cc.

266  {
267  std::sort(ohits.begin(),ohits.end(),sortOpHitByTime);
268 
269  std::vector<bool> isClust;
270  for (int i = 0; i < int(ohits.size()); i++) isClust.push_back(false);
271 
272  std::vector<int> neigh;
273 
274  for (int i = 0; i < int(ohits.size()); i++){
275 
276  if (isClust[i]) continue; // Don't base clusts off of clustered hits!
277  neigh.erase(neigh.begin(),neigh.end()); // Start from scratch every time
278 
279  // check nearby hits in time for coincidence
280  for (int j = i-1; j > 0; j--){
281  if ( Dist(ohits[i],ohits[j]) < fR0 && !isClust[j] ){
282  neigh.push_back(j);
283  }
284  if (abs(ohits[i]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2) break;
285  }
286  for (int j = i+1; j < int(ohits.size()); j++){
287  if ( Dist(ohits[i],ohits[j]) < fR0 && !isClust[j] ){
288  neigh.push_back(j);
289  }
290  if (abs(ohits[i]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2) break;
291  }
292  if (int(neigh.size())<fMinN) continue;
293  neigh.erase(neigh.begin(),neigh.end());
294 
295 
296  std::vector<int> cands;
297  for (int j = i; j < int(ohits.size()); j++){
298  if (isClust[j]) continue;
299  if (ohits[j]->PeakTimeAbs()-ohits[i]->PeakTimeAbs()>fBreakTime) break;
300  cands.push_back(j);
301  }
302  int centroidIdx = YZCentroid(ohits,cands);
303  if (centroidIdx < 0) // No centroid found
304  continue;
305 
306  art::Ptr<recob::OpHit> centroid = ohits[centroidIdx];
307  neigh.push_back(centroidIdx);
308 
309  std::vector<int> curN;
310  curN.push_back(centroidIdx);
311 
312  // check nearby hits in time for coincidence
313  for (int j = centroidIdx-1; j > 0; j--){
314  if ( Dist(ohits[centroidIdx],ohits[j]) < fR0 && !isClust[j] ){
315  neigh.push_back(j);
316  curN.push_back(j);
317  }
318  if (abs(ohits[centroidIdx]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2) break;
319  }
320  for (int j = centroidIdx+1; j < int(ohits.size()); j++){
321  if ( Dist(ohits[centroidIdx],ohits[j]) < fR0 && !isClust[j] ){
322  neigh.push_back(j);
323  curN.push_back(j);
324  }
325  if (abs(ohits[centroidIdx]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2) break;
326  }
327  double totPE = 0;
328  for (int idx : curN) totPE += ohits[idx]->PE();
329  if (int(curN.size())<fMinN) continue;
330 
331  // Loop through neighboring hits, chck if it's a core hit
332  while (neigh.size() > 0){
333  std::vector<int> curNeigh; curNeigh.push_back(neigh[0]);
334  for (int j = neigh[0]-1; j > 0; j--){
335  if ( Dist(ohits[neigh[0]],ohits[j]) < fR0 && !isClust[j] ){
336  curNeigh.push_back(j);
337  }
338  if (abs(ohits[centroidIdx]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2)
339  break;
340  }
341  for (int j = neigh[0]+1; j < int(ohits.size()); j++){
342  if ( Dist(ohits[neigh[0]],ohits[j]) < fR0 && !isClust[j] ){
343  curNeigh.push_back(j);
344  }
345  if (abs(ohits[centroidIdx]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2)
346  break;
347  }
348  // If this is a core point, add in all reachable hits to neighborhood
349  if (int(curNeigh.size())>=fMinN){
350  for (int cur : curNeigh){
351  if (std::find(curN.begin(),curN.end(),cur)==curN.end()){
352  curN.push_back(cur);
353  if (Dist(ohits[centroidIdx],ohits[cur]) < fR0)
354  neigh.push_back(cur);
355  }
356  }
357  }
358  neigh.erase(neigh.begin());
359  }
360 
361  if (int(curN.size())<fMinN) continue;
362 
363  // Time to make the OpFlash;
364  // Y-Z coordinates come from the centroid
365  int channelcentroid = centroid->OpChannel();
366  double xyzcentroid[3];
367  geo->OpDetGeoFromOpChannel(channelcentroid).GetCenter(xyzcentroid);
368  double yCenter = xyzcentroid[1];
369  double zCenter = xyzcentroid[2];
370  double tCenter = centroid->PeakTimeAbs();
371 
372 
373  // Now that we have centroid coordinates, include ana delayed light
374  for (int j = i; j < int(ohits.size()); j++){
375  if (std::find(curN.begin(),curN.end(),j)!=curN.end()) continue;
376  double r = YZDist(ohits[j],centroid);
377  if ( r < fRScale*fR0){
378  curN.push_back(j);
379  }
380  if (abs(ohits[j]->PeakTimeAbs()-tCenter)>fBreakTime) break;
381  }
382 
383  double finE = 0;
384  for (int idx : curN) finE += ohits[idx]->PE();
385 
386  // Grab the y-z information from the geometry
387  std::vector<double> ys;
388  std::vector<double> zs;
389  GetHitYZ(ohits,curN,ys,zs);
390 
391  double minT = std::numeric_limits<double>::max(); double maxT = -std::numeric_limits<double>::max();
392  double minY = 1e6; double maxY = -1e6;
393  double minZ = 1e6; double maxZ = -1e6;
394 
395  std::vector<double> PEs (geo->MaxOpChannel() + 1,0.0);
396  std::vector<double> PE2s (geo->MaxOpChannel() + 1,0.0);
397  double fastToTotal = 0;
398  for (int hIdx = 0; hIdx < int(ys.size()); hIdx++){
399  int cIdx = curN[hIdx];
400 
401  minT = std::min(minT,ohits[cIdx]->PeakTimeAbs());
402  maxT = std::max(maxT,ohits[cIdx]->PeakTimeAbs());
403  minY = std::min(minY,ys[hIdx]);
404  maxY = std::min(maxY,ys[hIdx]);
405  minZ = std::min(minZ,zs[hIdx]);
406  maxZ = std::min(maxZ,zs[hIdx]);
407  PEs[ohits[cIdx]->OpChannel()] += ohits[cIdx]->PE();
408  PE2s[ohits[cIdx]->OpChannel()] += ohits[cIdx]->PE();
409  fastToTotal += ohits[hIdx]->FastToTotal();
410  }
411  double yWidth = maxY-minY;
412  double zWidth = maxZ-minZ;
413 
414  double tot1 = 0;
415  double tot2 = 0;
416  for (double PE : PEs) tot1 += PE;
417  for (double PE : PE2s) tot2 += PE;
418 
419  // From OpFlashAlg
420  int Frame = ts.OpticalClock().Frame(tCenter - 18.1);
421  if (Frame == 0) Frame = 1;
422 
423  int BeamFrame = ts.OpticalClock().Frame(ts.TriggerTime());
424  bool InBeamFrame = false;
425  if (!(ts.TriggerTime() < 0)) InBeamFrame = (Frame == BeamFrame);
426 
427  double tWidth = (maxT-minT)/2;
428 
429  int OnBeamTime = 0;
430  if (InBeamFrame && (std::abs(tCenter) < fTrigCoinc)) OnBeamTime = 1;
431 
432  oflashes.emplace_back(tCenter,tWidth,tCenter,Frame,
433  PEs,InBeamFrame,OnBeamTime,fastToTotal,
434  yCenter,yWidth,zCenter,zWidth);
435  assoc.emplace_back(curN);
436 
437 
438  // And finally, indicate that current hits have been clustered
439  for (int cur : curN) isClust[cur] = true;
440 
441  }
442  }
void GetHitYZ(std::vector< art::Ptr< recob::OpHit > >, std::vector< int >, std::vector< double > &, std::vector< double > &) const
int YZCentroid(std::vector< art::Ptr< recob::OpHit > >, std::vector< int >) const
double YZDist(art::Ptr< recob::OpHit >, art::Ptr< recob::OpHit >) const
T abs(T value)
bool sortOpHitByTime(const art::Ptr< recob::OpHit > &left, const art::Ptr< recob::OpHit > &right)
static int max(int a, int b)
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
double PeakTimeAbs() const
Definition: OpHit.h:65
static constexpr double zs
Definition: Units.h:102
double Dist(art::Ptr< recob::OpHit >, art::Ptr< recob::OpHit >) const
static constexpr double ys
Definition: Units.h:103
LArSoft geometry interface.
Definition: ChannelGeo.h:16
int OpChannel() const
Definition: OpHit.h:62
double opdet::OpSlicer::Dist ( art::Ptr< recob::OpHit a,
art::Ptr< recob::OpHit b 
) const
private

Definition at line 218 of file OpSlicer_module.cc.

220  {
221  double r = YZDist(a,b);
222  return sqrt(pow((a->PeakTime()-b->PeakTime())/fTScale,2)+pow(r/fRScale,2));
223  }
constexpr T pow(T x)
Definition: pow.h:72
double PeakTime() const
Definition: OpHit.h:64
double YZDist(art::Ptr< recob::OpHit >, art::Ptr< recob::OpHit >) const
void opdet::OpSlicer::endJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 145 of file OpSlicer_module.cc.

146  {
147  }
void opdet::OpSlicer::GetHitYZ ( std::vector< art::Ptr< recob::OpHit > >  ohits,
std::vector< int >  curN,
std::vector< double > &  ys,
std::vector< double > &  zs 
) const
private

Definition at line 246 of file OpSlicer_module.cc.

250  {
251  for (int cur : curN){
252  art::Ptr<recob::OpHit> oh = ohits[cur];
253  int channel = oh->OpChannel();
254  double xyz[3];
255  geo->OpDetGeoFromOpChannel(channel).GetCenter(xyz);;
256  ys.push_back(xyz[1]);
257  zs.push_back(xyz[2]);
258  }
259  }
uint8_t channel
Definition: CRTFragment.hh:201
LArSoft geometry interface.
Definition: ChannelGeo.h:16
int OpChannel() const
Definition: OpHit.h:62
void opdet::OpSlicer::produce ( art::Event evt)
virtual

Implements art::EDProducer.

Definition at line 150 of file OpSlicer_module.cc.

151  {
152 
153  // These are the storage pointers we will put in the event
154  std::unique_ptr< std::vector< recob::OpFlash > >
155  flashPtr(new std::vector< recob::OpFlash >);
156  std::unique_ptr< art::Assns< recob::OpFlash, recob::OpHit > >
158 
159  // This will keep track of what flashes will assoc to what ophits
160  // at the end of processing
161  std::vector< std::vector< int > > assocList;
162 
163  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(evt);
164 
165  // Get OpHits from the event
166  auto opHitHandle = evt.getHandle< std::vector< recob::OpHit > >(fOpHitModuleLabel);
167 
168  std::vector< art::Ptr<recob::OpHit> > ohits;
169  for (int i = 0; i < int(opHitHandle->size()); i++){
170  art::Ptr<recob::OpHit> opHitPtr(opHitHandle,i);
171  ohits.push_back(opHitPtr);
172  }
173 
174 
175  // Run the clustering
176  ClusterHits(ohits, *flashPtr, assocList, clockData);
177 
178  // Make the associations which we noted we need
179  for (size_t i = 0; i != assocList.size(); ++i)
180  {
181  art::PtrVector< recob::OpHit > opHitPtrVector;
182  for (int const& hitIndex : assocList.at(i))
183  {
184  art::Ptr< recob::OpHit > opHitPtr(opHitHandle, hitIndex);
185  opHitPtrVector.push_back(opHitPtr);
186  }
187 
188  util::CreateAssn(*this, evt, *flashPtr, opHitPtrVector,
189  *(assnPtr.get()), i);
190  }
191 
192  // Store results into the event
193  evt.put(std::move(flashPtr));
194  evt.put(std::move(assnPtr));
195  }
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
def move(depos, offset)
Definition: depos.py:107
void ClusterHits(std::vector< art::Ptr< recob::OpHit >>, std::vector< recob::OpFlash > &, std::vector< std::vector< int > > &, detinfo::DetectorClocksData const &) const
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
std::string fOpHitModuleLabel
void opdet::OpSlicer::reconfigure ( fhicl::ParameterSet const &  pset)

Definition at line 118 of file OpSlicer_module.cc.

119  {
120 
121  // Indicate that the Input Module comes from .fcl
122  fOpHitModuleLabel = pset.get< std::string >("OpHitModuleLabel");
123 
124  fTScale = pset.get<double>("TScale");
125  fRScale = pset.get<double>("RScale");
126  fR0 = pset.get<double>("R0");
127  fBreakTime = pset.get<double>("BreakTime");
128  fMinN = pset.get<int> ("MinN");
129 
130  fTrigCoinc = pset.get< double >("TrigCoinc");
131 
132  }
std::string string
Definition: nybbler.cc:12
std::string fOpHitModuleLabel
int opdet::OpSlicer::YZCentroid ( std::vector< art::Ptr< recob::OpHit > >  ohits,
std::vector< int >  curN 
) const
private

Definition at line 225 of file OpSlicer_module.cc.

227  {
228  double maxDens = 0;
229  int maxIdx = -1;
230  for (int i = 0; i < int(curN.size()); i++){
231  double dens = 0;
232  for (int j = 0; j < int(curN.size()); j++){
233  double r = YZDist(ohits[i],ohits[j]);
234  r = sqrt(r)/100;
235  dens += ohits[curN[j]]->PE() * exp(-r*r);
236  }
237  if (dens > maxDens){
238  maxDens = dens; maxIdx = curN[i];
239  }
240  }
241  return maxIdx;
242  }
double YZDist(art::Ptr< recob::OpHit >, art::Ptr< recob::OpHit >) const
double opdet::OpSlicer::YZDist ( art::Ptr< recob::OpHit a,
art::Ptr< recob::OpHit b 
) const
private

Definition at line 197 of file OpSlicer_module.cc.

199  {
200  // First need to ask the geometry for the y-z location of both OpHits
201  int channela = a->OpChannel();
202  double xyza[3];
203  geo->OpDetGeoFromOpChannel(channela).GetCenter(xyza);;
204  double ay = xyza[1];
205  double az = xyza[2];
206 
207  int channelb = b->OpChannel();
208  double xyzb[3];
209  geo->OpDetGeoFromOpChannel(channelb).GetCenter(xyzb);;
210  double by = xyzb[1];
211  double bz = xyzb[2];
212 
213  double r2 = pow((ay-by),2) + pow((az-bz),2);
214 
215  double r = sqrt(r2);
216  return r;
217  }
constexpr T pow(T x)
Definition: pow.h:72
LArSoft geometry interface.
Definition: ChannelGeo.h:16
int OpChannel() const
Definition: OpHit.h:62

Member Data Documentation

double opdet::OpSlicer::fBreakTime
private

Definition at line 87 of file OpSlicer_module.cc.

int opdet::OpSlicer::fMinN
private

Definition at line 88 of file OpSlicer_module.cc.

std::string opdet::OpSlicer::fOpHitModuleLabel
private

Definition at line 82 of file OpSlicer_module.cc.

double opdet::OpSlicer::fR0
private

Definition at line 86 of file OpSlicer_module.cc.

double opdet::OpSlicer::fRScale
private

Definition at line 85 of file OpSlicer_module.cc.

double opdet::OpSlicer::fTrigCoinc
private

Definition at line 90 of file OpSlicer_module.cc.

double opdet::OpSlicer::fTScale
private

Definition at line 84 of file OpSlicer_module.cc.

art::ServiceHandle<geo::Geometry> opdet::OpSlicer::geo
private

Definition at line 92 of file OpSlicer_module.cc.


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