Recluster_module.cc
Go to the documentation of this file.
1 
2 #include <fstream>
3 #include <cstdlib>
4 #include <vector>
5 #include <string>
6 #include "TVector3.h"
7 
8 // ART includes
18 
19 // NOvA includes
20 #include "Geometry/Geometry.h"
22 #include "RecoBase/CellHit.h"
23 #include "RecoBase/RecoHit.h"
24 #include "RecoBase/Cluster.h"
25 #include "RecoBase/Shower.h"
26 #include "RecoBase/FilterList.h"
27 #include "Utilities/AssociationUtil.h"
29 #include "RecoBase/Vertex.h"
30 
31 namespace slid{
32 
33  class Recluster : public art::EDProducer
34  {
35 
36  public:
37 
38  Recluster( fhicl::ParameterSet const& pset );
39  ~Recluster();
40 
41  void produce(art::Event& evt);
42  void reconfigure(const fhicl::ParameterSet& pset);
43  void beginRun(art::Run& run);
44 
45  private:
46 
47  std::string fSliceLabel; ///< Label of slicer module
48  bool fProngInput; ///< Lable of prong producing module
49  std::string fProngLabel; ///< True = Prong input, False = Track input
50  std::string fVertexLabel; ///< Lable of vertex producing module
51  std::string fProng3DInstLabel; ///< Label of instance label of 3d prongs
52  std::string fProng2DInstLabel; ///< Label of instance label of 2d prongs
53  std::string fTrack3DLabel; ///< Label of 3d tracks
54  std::string fTrack2DLabel; ///< Label of 2d tracks
55  std::vector<std::string> fFilterLabels; ///< Labels of filter lists to obey
56 
57  bool fRespectFilter; ///< Skip slices that failed filter?
58  bool fRecluster; ///< Should shower reclustering be performed?
59  bool fDeconvolve; ///< Should cell-energy deconvolution be performed?
60  bool fConsider2D; ///< Should 2-D prongs be considered for reclustering/deconvolution?
61  bool fApplyPEThreshold; ///< If this parameter is true, cellhits with PE below fPEThreshold will not be reclustered
62  float fMaxPlane; ///< The extent of the reclustering cylinder, beyond prong stop point
63  float fMIPPlane; ///< The extent of the MIP cylinder from the prong start point
64  float fMaxRadius; ///< The radius of the reclustering cylinder, in cell widths
65  float fMIPRadius; ///< The radius of the MIP cylinder, in cell widths
66 
67  float fPEThreshold; ///< If a cell PE is less than this value, it will not be added to the shower.
68 
69 
70  DeconvolveAlg* fDeconvolveAlg; ///< An object of DeconvolveAlg
71  fhicl::ParameterSet fDeconvolveAlgPSet; ///< ParameterSet of DecolvolveAlg
72  art::ServiceHandle<geo::Geometry> fGeom; ///< A private handle to geometry service
73 
74  };// end Recluster class definition
75 
76 
77  /**********************************************************************/
78  ////////////////////////////////////////////////////////////////////////
79 
80 
82  : fDeconvolveAlg(0)
83  , fDeconvolveAlgPSet(pset.get< fhicl::ParameterSet >("DeconvolveAlgPSet") )
84  {
85  this->reconfigure(pset);
86 
87  // Produces a vector of reclustered showers and their
88  // association with slices
89  produces< std::vector <rb::Shower> >();
90  produces< art::Assns <rb::Shower, rb::Cluster> >();
91  produces< art::Assns <rb::Shower, rb::Prong> >();
92 
93  }// End of constructor
94 
95 
96  /**********************************************************************/
97  ////////////////////////////////////////////////////////////////////////
98 
100  {
101  if( fDeconvolveAlg )
102  delete fDeconvolveAlg;
103  }
104 
105  /**********************************************************************/
106  ////////////////////////////////////////////////////////////////////////
107 
109  {
110  fSliceLabel = pset.get< std::string >("SliceLabel");
111  fProngInput = pset.get< bool >("ProngInput");
112  fProngLabel = pset.get< std::string >("ProngLabel");
113  fVertexLabel = pset.get< std::string >("VertexLabel");
114  fProng3DInstLabel = pset.get< std::string >("Prong3DInstLabel");
115  fProng2DInstLabel = pset.get< std::string >("Prong2DInstLabel");
116  fTrack3DLabel = pset.get< std::string >("Track3DLabel");
117  fTrack2DLabel = pset.get< std::string >("Track2DLabel");
118 
119 
120  fMaxPlane = pset.get< float >("MaxPlane");
121  fMIPPlane = pset.get< float >("MIPPlane");
122  fMaxRadius = pset.get< float >("MaxRadius");
123  fMIPRadius = pset.get< float >("MIPRadius");
124  fPEThreshold = pset.get< float >("PEThreshold");
125 
126  fRespectFilter = pset.get< bool >("RespectFilter");
127  fRecluster = pset.get< bool >("Recluster");
128  fDeconvolve = pset.get< bool >("Deconvolve");
129  fConsider2D = pset.get< bool >("Consider2D");
130  fApplyPEThreshold = pset.get< bool >("ApplyPEThreshold");
131  fFilterLabels = pset.get< std::vector<std::string> >("FilterLabels");
132  }
133 
134  /*********************************************************************/
135  //////////////////////////////////////////////////////////////////////
136 
138  {
140  }
141 
142  /**********************************************************************/
143  ////////////////////////////////////////////////////////////////////////
144 
146  {
147 
148  std::unique_ptr< std::vector< rb::Shower > >
149  reclustCol( new std::vector< rb::Shower > );
150  std::unique_ptr< art::Assns< rb::Shower, rb::Cluster > >
151  sliceAssassin( new art::Assns< rb::Shower, rb::Cluster > );
152  // TODO: This is strictly temporary. Once validation
153  // is over, we won't need prong-shower association any
154  // longer and shouldn't be produced.
155  std::unique_ptr< art::Assns< rb::Shower, rb::Prong > >
156  prongAssassin( new art::Assns< rb::Shower, rb::Prong > );
157 
158  // grab the slices in this event
160  evt.getByLabel( fSliceLabel, slices );
161 
162  std::vector< rb::Shower > showerVec;
163  // grab the prong-slice associations in this event
165  fmprong3D(slices, evt, art::InputTag(fProngLabel, fProng3DInstLabel));
167  fmprong2D(slices, evt, art::InputTag(fProngLabel, fProng2DInstLabel));
168 
170  fmtrack3D(slices, evt, fTrack3DLabel);
172  fmtrack2D(slices, evt, fTrack2DLabel);
173 
174 
175 
176  // grab the vertex-slice associations in this event
177  art::FindManyP<rb::Vertex> fmvertex(slices, evt, fVertexLabel);
178 
179  // If there are no vertices or no prongs in this event
180  // move on.
181 
182  if( ((fProngInput==true&&fmprong3D.isValid())||(fProngInput==false&&fmtrack3D.isValid())) || fmvertex.isValid() ) {
183 
184  int nSli = slices->size();
185  for(int iSli = 0; iSli < nSli; ++iSli){
186  //if there was nue preselection, skip slice if it failed
187  if(rb::IsFiltered(evt,slices,iSli, fFilterLabels)
188  && fRespectFilter )
189  continue;
190 
191  art::Ptr<rb::Cluster> thisSlice(slices, iSli);
192 
193  // if this slice is noise skip it
194  if( thisSlice->IsNoise())
195  continue;
196 
197 // std::vector<art::Ptr<rb::Prong> > prongs = fmprong3D.at(iSli);
198 
199  std::vector<art::Ptr<rb::Track> > tracks3D;
200  tracks3D.clear();
201  if(fProngInput==false) tracks3D = fmtrack3D.at(iSli);
202 
203 
204  std::vector<art::Ptr<rb::Prong> > prongs;
205  prongs.clear();
206  if(fProngInput==true){//prong input
207  std::vector<art::Ptr<rb::Prong> > prongs3D = fmprong3D.at(iSli);
208  for (unsigned int ip=0; ip<prongs3D.size(); ++ip){
209  prongs.push_back(prongs3D[ip]);
210  }
211  //if instance labels are in use, combine 2d and 3d
212  if ((fProng3DInstLabel != "")&&(fProng2DInstLabel != "")&&fmprong2D.isValid()){
213  std::vector<art::Ptr<rb::Prong> > prongs2D = fmprong2D.at(iSli);
214  for (unsigned int ip=0; ip<prongs2D.size(); ++ip){
215  prongs.push_back(prongs2D[ip]);
216  }
217  }
218  }
219 
220  if(fProngInput==false){// track input, convert tracks into prongs
221  std::vector<art::Ptr<rb::Track> > tracks3D = fmtrack3D.at(iSli);
222  for (unsigned int ip=0; ip<tracks3D.size(); ++ip){
223  prongs.push_back(tracks3D[ip]);
224  }
225  if(fmtrack2D.isValid()){
226  std::vector<art::Ptr<rb::Track> > tracks2D = fmtrack2D.at(iSli);
227  for (unsigned int ip=0; ip<tracks2D.size(); ++ip){
228  prongs.push_back(tracks2D[ip]);
229  }
230  }
231  }
232 
233  int nPro = prongs.size();
234 
235  std::vector< rb::Prong > newProngs;
236  std::vector< int > assProng;
237 
238  std::vector<art::Ptr<rb::Vertex> > vertices = fmvertex.at(iSli);
239  if( vertices.empty() )
240  continue;
241  // Right now we only use fuzzyk for vertices which makes
242  // one vertex per slice. If there comes along another
243  // vertex finder that makes more than one per slice,
244  // the hope is that the 0th one will be the primary vertex.
245  art::Ptr<rb::Vertex> vtx = vertices[0];
246 
247  // Make a list of hits in slice that were not in
248  // any prong. These will be subject to a PE cut
249  // since hits not on prongs are likely to be noise
250 
251  rb::Cluster orphans = *(thisSlice);
252  if( fApplyPEThreshold ){
253  for( int iPro = 0; iPro < nPro; ++iPro){
254  rb::Prong prong = *(prongs[iPro]);
255  if( prong.Is2D() && !fConsider2D)
256  continue;
257  int nc = prong.NCell();
258  for( int ic = 0; ic < nc; ++ic)
259  orphans.RemoveHit( prong.Cell(ic));
260  }
261 
262  for(int io = 0; io < (int)orphans.NCell(); ++io)
263  if( orphans.Cell(io)->PE() > fPEThreshold ){
264  orphans.RemoveHit( io );
265  io = io -1;
266  }
267  }
268  // We should now have a list of below threshold orphan hits
269 
270  for( int iPro = 0; iPro < nPro; ++iPro){
271 
272  rb::Prong thisProng = *(prongs[iPro]);
273  if( thisProng.Is2D() && !fConsider2D )
274  continue;
275 
276  if( fRecluster ){
277 
278  const rb::Cluster showerClust = (rb::Cluster)thisProng;
279 
280  // Get the cluster that is the slice minus this prong
281  rb::Cluster lesserSlice = (showerClust.NCell() < thisSlice->NCell()) ? thisSlice->Exclude( &showerClust ) : rb::Cluster();
282  if( fApplyPEThreshold ){
283  int norphans = orphans.NCell();
284  for(int io = 0; io < norphans; ++io)
285  lesserSlice.RemoveHit( orphans.Cell(io) );
286  }
287 
288  // Get all the hits in this cluster and sort by plane-
289  art::PtrVector<rb::CellHit> sliceHits = lesserSlice.AllCells();
290 
291  rb::SortByPlaneAndCell( sliceHits );
292 
293  // Reclustering is done to grab electron hits that may
294  // have been missed by the prong making algorithm.
295  // To recluster, a series of two collinear cylinders is
296  // defined. The first cylinder captures the MIP behav-
297  // -iour of the electron; the second cylinder captures the
298  // EM shower. Their radii and lengths are set in .fcl file
299 
300 
301  TVector3 dir = thisProng.Dir();
302  TVector3 start = thisProng.Start();
303  TVector3 stop = start + thisProng.TotalLength()*dir;
304  int stopPlane, startPlane, vtxPlane;
305 
306  if( (start - vtx->GetXYZ()).Mag() > (stop - vtx->GetXYZ()).Mag() ){
307  // vertex is closer to the stop point. Reverse
308  // the shower direction.
309  TVector3 temp = start;
310  start = stop;
311  stop = temp;
312  thisProng.SetDir( -1*dir);
313  thisProng.SetStart( start );
314  }
315 
316  startPlane = thisProng.MinPlane();
317  stopPlane = thisProng.MaxPlane();
318 
319 
320  // There is a better way of determining shower start and stop planes
321  // than just using the cluster min and max planes.
322  // And that is to use geo::CellID as follows. However, the following
323  // is being kept commented, to agree with RecoJMShower for now. (2014/07/11)
324 
325  // const geo::CellUniqueId cida = fGeom->CellId( start.X(), start.Y(),
326  // start.Z(), 1, 1, 0, 1 );
327  // if( cida )
328  // startPlane = fGeom->getPlaneID( cida );
329  // else
330  // startPlane = thisProng.MinPlane();
331 
332  // const geo::CellUniqueId cidb = fGeom->CellId( stop.X(), stop.Y(),
333  // stop.Z(), 1, 1, 0, 1 );
334  // if( cidb )
335  // stopPlane = fGeom->getPlaneID( cidb );
336  // else
337  // stopPlane = thisProng.MaxPlane();
338 
339  const geo::CellUniqueId cidv = fGeom->CellId( vtx->GetX(),
340  vtx->GetY(),
341  vtx->GetZ(),
342  1, 1, 0, 1 );
343 
344  if( cidv )
345  vtxPlane = fGeom->getPlaneID( cidv );
346  else
347  vtxPlane = thisProng.MinPlane();
348 
349  if( abs(vtxPlane - startPlane) > abs(vtxPlane-stopPlane))
350  std::swap(startPlane, stopPlane);
351 
352 
353  float cellW = 2*fGeom->Plane(1)->Cell(1)->HalfW();
354 
355  for( auto hit : sliceHits ){
356 
357  int thisPlane = hit->Plane();
358  int diffFromStart = thisPlane - startPlane;
359 
360  if( startPlane > stopPlane ){
361  diffFromStart = -diffFromStart;
362  }
363 
364  if( thisPlane < std::min(startPlane, stopPlane) ||
365  thisPlane > std::max(startPlane, stopPlane) )
366  continue;
367 
368  float distToCore = fDeconvolveAlg->DistanceToCore( hit,
369  thisProng );
370 
371  if( distToCore == 9999 )
372  continue;
373 
374  if( (diffFromStart < fMIPPlane &&
375  distToCore <= fMIPRadius*cellW) ||
376  (diffFromStart >= fMIPPlane &&
377  distToCore <= fMaxRadius*cellW))
378  thisProng.Add( hit );
379 
380  }// end loop over relevant cell hits
381  }// end if recluster
382  newProngs.push_back( thisProng );
383  assProng.push_back( iPro );
384  }// end loop over prongs
385 
386  std::vector< float > totpe;
387  for( int iPro = 0; iPro < nPro; ++iPro){
388  totpe.push_back( prongs[iPro]->TotalPE() );
389  }
390  if( fDeconvolve )
391  fDeconvolveAlg->GetWeights( newProngs, totpe );
392 
393  int nNewShower = newProngs.size();
394 
395  for( int i = 0; i < nNewShower; ++i){
396  rb::Shower newShower(newProngs[i]);
397  newShower.SetID(i + (100*iSli) );
398  newShower.SetTotalLength( newProngs[i].TotalLength() );
399  reclustCol->push_back( newShower );
400  util::CreateAssn( *this, evt, *reclustCol,
401  thisSlice ,*sliceAssassin);
402  util::CreateAssn( *this, evt, *reclustCol,
403  prongs[ assProng[i] ] ,*prongAssassin);
404  }
405 
406  }// end loop over slices
407 
408  }// end if slice-vertex or slice-prong associations are invalid
409 
410 
411  evt.put( std::move(reclustCol));
412  evt.put( std::move(sliceAssassin));
413  evt.put( std::move(prongAssassin));
414  }// End of producer
415 
416 
417  /**********************************************************************/
418  ////////////////////////////////////////////////////////////////////////
419 
421 
422 }// end of namespace slid
T max(const caf::Proxy< T > &a, T b)
bool fConsider2D
Should 2-D prongs be considered for reclustering/deconvolution?
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
art::ServiceHandle< geo::Geometry > fGeom
A private handle to geometry service.
bool fRespectFilter
Skip slices that failed filter?
bool fDeconvolve
Should cell-energy deconvolution be performed?
Preforms energy deconvolution for cells that may belong to multple clusters. It does so...
bool fApplyPEThreshold
If this parameter is true, cellhits with PE below fPEThreshold will not be reclustered.
bool Is2D() const
Definition: Cluster.h:96
float DistanceToCore(art::Ptr< rb::CellHit > cHit, rb::Prong &prong)
double GetX() const
Definition: Vertex.h:23
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
std::string fSliceLabel
Label of slicer module.
void beginRun(art::Run &run)
double HalfW() const
Definition: CellGeo.cxx:191
A collection of associated CellHits.
Definition: Cluster.h:47
void abs(TH1 *hist)
void RemoveHit(const art::Ptr< rb::CellHit > hit)
Remove hit from current cluster.
Definition: Cluster.cxx:290
TString ip
Definition: loadincs.C:5
double GetY() const
Definition: Vertex.h:24
const PlaneGeo * Plane(unsigned int i) const
rb::Cluster Exclude(const rb::Cluster *excl) const
Create a cluster from this one, but with the hits of excl removed.
Definition: Cluster.cxx:233
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
void produce(art::Event &evt)
void SetTotalLength(double length)
Set length of a shower.
Definition: Shower.cxx:49
virtual void SetStart(TVector3 start)
Definition: Prong.cxx:75
fhicl::ParameterSet fDeconvolveAlgPSet
ParameterSet of DecolvolveAlg.
float fPEThreshold
If a cell PE is less than this value, it will not be added to the shower.
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
double GetZ() const
Definition: Vertex.h:25
Definition: Run.h:31
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
bool fRecluster
Should shower reclustering be performed?
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
std::string fTrack2DLabel
Label of 2d tracks.
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Definition: Prong.cxx:186
TVector3 GetXYZ() const
Definition: Vertex.cxx:45
virtual void SetDir(TVector3 dir)
Definition: Prong.cxx:104
T get(std::string const &key) const
Definition: ParameterSet.h:231
int getPlaneID(const CellUniqueId &id) const
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
int evt
float PE() const
Definition: CellHit.h:42
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
Recluster(fhicl::ParameterSet const &pset)
void SetID(int id)
Definition: Cluster.h:74
std::string fProngLabel
True = Prong input, False = Track input.
Vertex location in position and time.
float fMaxPlane
The extent of the reclustering cylinder, beyond prong stop point.
Definition: run.py:1
void reconfigure(const fhicl::ParameterSet &pset)
void GetWeights(std::vector< rb::Prong > &prong, const std::vector< float > &totpe)
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
std::string fTrack3DLabel
Label of 3d tracks.
A Cluster with defined start position and direction.
Definition: Prong.h:19
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::string fProng2DInstLabel
Label of instance label of 2d prongs.
Definition: structs.h:12
A rb::Prong with a length.
Definition: Shower.h:18
float fMIPPlane
The extent of the MIP cylinder from the prong start point.
float Mag() const
std::string fProng3DInstLabel
Label of instance label of 3d prongs.
float fMaxRadius
The radius of the reclustering cylinder, in cell widths.
enum BeamMode nc
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
Build slid::LID objects to store electron ID, if asked for, otherwise, calculate LID info and make av...
Definition: FillPIDs.h:13
float fMIPRadius
The radius of the MIP cylinder, in cell widths.
std::string fVertexLabel
Lable of vertex producing module.
T min(const caf::Proxy< T > &a, T b)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
std::vector< std::string > fFilterLabels
Labels of filter lists to obey.
Encapsulate the geometry of one entire detector (near, far, ndos)
DeconvolveAlg * fDeconvolveAlg
An object of DeconvolveAlg.
void SortByPlaneAndCell(std::vector< art::Ptr< rb::CellHit > > &c)
Definition: CellHit.cxx:155
bool fProngInput
Lable of prong producing module.
enum BeamMode string