FindParent_module.cc
Go to the documentation of this file.
1 // C/C++ & root includes
2 #include <iostream>
3 #include "vector"
4 #include <TTree.h>
5 #include "TH1D.h"
6 #include <map>
7 
8 // ART includes
9 
23 
24 //NOvA includes
25 #include "Utilities/AssociationUtil.h"
29 #include "ReMId/ReMId.h"
30 #include "NumuEnergy/NumuE.h"
31 #include "RecoBase/Cluster.h"
32 #include "RecoBase/Track.h"
33 #include "MCCheater/BackTracker.h"
34 #include "SummaryData/POTSum.h"
35 
36 namespace murem
37 {
38 
39  namespace
40  {
41  struct FindParentParams
42  {
43  template<class T> using Atom = fhicl::Atom<T>;
44  template<class T> using Table = fhicl::Table<T>;
45  using Comment = fhicl::Comment;
46  using Name = fhicl::Name;
47 
48  Atom<std::string> mrccLabel
49  {
50  Name("MrccLabel"),
51  Comment("Label of the mrcc raw digit module")
52  };
53  Atom<std::string> mrClusterLabel
54  {
55  Name("MRClusterLabel"),
56  Comment("Label of mrcc slicer module")
57  };
58  Atom<std::string> mrHitLabel
59  {
60  Name("MRHitLabel"),
61  Comment("Label of the mrcc calhit module")
62  };
63  Atom<std::string> hitLabel
64  {
65  Name("HitLabel"),
66  Comment("Label of the parent calhit module")
67  };
68  Atom<std::string> clusterLabel
69  {
70  Name("ClusterLabel"),
71  Comment("Label of the parent slicer module")
72  };
73 
74  };// end of parameters' struct
75 
76  }
77 
78  class FindParent : public art::EDProducer
79  {
80  public:
81 
83  explicit FindParent(const Parameters& params);
84  virtual ~FindParent();
85 
86  void produce(art::Event& evt);
87 
88  private:
89 
90  int MatchToOrigSlice(art::Ptr<rb::Cluster> mrSlice,
91  art::Handle< std::vector<rb::Cluster> > origSlices,
92  std::vector< art::Ptr< rb::CellHit> > mrHits,
93  float& maxPur,
94  float& maxEff);
95 
96  FindParentParams fParams;
97  };
98 
99  /*********************************************************************************************/
100  ////////////////////////////////////////////////////////////////////////
101 
103  : fParams(params())
104  {
105 
106  produces<std::vector<murem::MRCCParent> >();
107  produces< art::Assns<murem::MRCCParent, rb::Cluster> >();
108  }
109 
110  /*********************************************************************************************/
111  ////////////////////////////////////////////////////////////////////////
112 
113 
115  {
116  }
117 
118 
119  /*********************************************************************************************/
120  ////////////////////////////////////////////////////////////////////////
122  {
123 
124  std::unique_ptr< std::vector< murem::MRCCParent > >
125  mrccParents(new std::vector< murem::MRCCParent >);
126 
127  std::unique_ptr< art::Assns< murem::MRCCParent, rb::Cluster > >
129 
130  std::vector< murem::MRCCParent> mrccParentsVec;
131 
133  evt.getByLabel(fParams.clusterLabel(), slicelist);
134  if (slicelist->empty()) {
135  mf::LogWarning("FindParent") << "Error retrieving slice list" ;
136  return;
137  }
138 
139 
141  evt.getByLabel(fParams.mrClusterLabel(), mrslicelist);
142  if (mrslicelist->empty()) {
143  mf::LogWarning("FindParent") << "Error retrieving mrcc slice list" ;
144  return;
145  }
146 
148  evt.getByLabel(fParams.mrccLabel(), removedMuonCol);
149 
151  for(int i = 0; i < (int)removedMuonCol->size(); i++){
152  art::Ptr<simb::MCParticle> mu( removedMuonCol, i );
153  removedMuons.push_back( mu );
154  }
155 
156 
157  // Grab all the MRCC cellhits. Need this for parent slice efficiency calculation
158 
160  evt.getByLabel(fParams.mrHitLabel(), mrHitlist);
161 
162  std::vector< art::Ptr<rb::CellHit> > mrAllhits;
163  for(size_t h = 0; h < mrHitlist->size(); ++h){
164  art::Ptr<rb::CellHit> hit(mrHitlist, h);
165  mrAllhits.push_back(hit);
166  }
167 
168  int nSlicesMrcc = mrslicelist->size();
169 
170  for(int mrSli = 0; mrSli < nSlicesMrcc; mrSli++ ){
171 
172  art::Ptr<rb::Cluster> thisMrSlice(mrslicelist, mrSli);
173  if( thisMrSlice->IsNoise() )
174  continue;
175 
176  float slicePur = -5;
177  float sliceEff = -5;
178 
179 
180  int parentSliceIndex = MatchToOrigSlice( thisMrSlice, slicelist, mrAllhits,
181  slicePur, sliceEff);
182  if (parentSliceIndex == -5 )
183  continue;
184 
185  MRCCParent mrccParent;
186  mrccParent.SetParentPur( slicePur );
187  mrccParent.SetParentEff( sliceEff );
188  mrccParent.SetParentIndex( parentSliceIndex );
189 
190  mrccParents->push_back( mrccParent );
191  util::CreateAssn( *this, evt, *mrccParents, thisMrSlice ,*sliceAssassin);
192 
193 
194  }// end loop over mrcc slices
195 
196 
197  evt.put( std::move(mrccParents) );
198  evt.put( std::move(sliceAssassin) );
199  }// End of produce method
200 
201 
202 
203 
204 
205  /*********************************************************************************************/
206  ////////////////////////////////////////////////////////////////////////
207 
209  art::Handle<std::vector<rb::Cluster> > origSlices,
210  std::vector< art::Ptr<rb::CellHit> > mrHits,
211  float& maxPur,
212  float& maxEff)
213  {
214 
215  SortByPlane(mrHits);
216 
217  auto minPlane = [](art::Ptr<rb::CellHit> hit, int plane ){ return (hit->Plane() < plane); };
218  auto maxPlane = [](int plane, art::Ptr<rb::CellHit> hit ){ return (hit->Plane() >= plane); };
219 
220  int nSlices = origSlices->size();
221  int nMrCells = mrccSlice->NCell();
222  int mrMin = mrccSlice->MinPlane();
223  int mrMax = mrccSlice->MaxPlane();
224  int matchedSliceIndex = -5;
225  maxPur = -1;
226 
227  for( int sli = 0; sli < nSlices; sli++ ){
228 
229  art::Ptr<rb::Cluster> thisSlice( origSlices, sli);
230 
231  // Don't match mrcc slices to original noise slice, even
232  // if it legitimate. Noise slice info is not saved in
233  // cafs and therefore, causes trouble at the time of analysis
234  if( thisSlice->IsNoise() )
235  continue;
236 
237  int origMin = thisSlice->MinPlane();
238  int origMax = thisSlice->MaxPlane();
239 
240 
241  if( !( ( mrMin >= origMin && mrMin <= origMax ) ||
242  ( mrMax >= origMin && mrMax <= origMax ) ) )
243  continue;
244 
245 
246  int nCells = thisSlice->NCell();
247 
248  int cellCount = 0;
249  for(int h = 0; h < nMrCells; h++){
250  art::Ptr<rb::CellHit> mhit = mrccSlice->Cell(h);
251  int plane = mhit->Plane();
252  int cell = mhit->Cell();
253  int tdc = mhit->TDC();
254 
255  for(int oh = 0; oh < nCells; oh++){
256  art::Ptr<rb::CellHit> hit = thisSlice->Cell(oh);
257  if( plane == hit->Plane() &&
258  cell == hit->Cell() &&
259  tdc == hit->TDC())
260  cellCount++;
261  }// end loop over orig slice hits
262 
263  }// end loop over hits in mrccSlice
264 
265 
266  float pur = ((float)cellCount)/nMrCells;
267 
268  if( pur > maxPur ){
269  matchedSliceIndex = sli;
270  maxPur = pur;
271  }
272 
273  }// end loop over the orig slices
274 
275  // Once the original slice for which the mrcc slice has the max
276  // purity is found, determine the efficiency of mrcc slice wrt
277  // hits in the original matched slice minus the muon hits.
278 
279  if(matchedSliceIndex == -5)
280  return matchedSliceIndex;
281 
282 
283  art::Ptr<rb::Cluster> origSlice( origSlices, matchedSliceIndex );
284  rb::Cluster nonMuonClust = *(origSlice);
285  int nCells = origSlice->NCell();
286  int minP = origSlice->MinPlane();
287  int maxP = origSlice->MaxPlane();
288 
289 
290  art::PtrVector<rb::CellHit>::iterator mrhitIter, mrhitBegin, mrhitEnd;
291  mrhitEnd = std::upper_bound( mrHits.begin(), mrHits.end(),
292  maxP, maxPlane );
293  mrhitBegin = std::lower_bound( mrHits.begin(), mrHits.end(),
294  minP, minPlane );
295 
296 
297  for(int oh = 0; oh < nCells; oh++){
298  const art::Ptr<rb::CellHit> hit = origSlice->Cell(oh);
299  int plane = hit->Plane();
300  int cell = hit->Cell();
301  int tdc = hit->TDC();
302  bool matchedHit = false;
303 
304  for( mrhitIter = mrhitBegin; mrhitIter != mrhitEnd; ++mrhitIter){
305  const art::Ptr<rb::CellHit> mrHit = *mrhitIter;
306  if( mrHit->Plane() == plane &&
307  mrHit->Cell() == cell &&
308  mrHit->TDC() == tdc )
309  matchedHit = true;
310 
311  }// end loop over mrcc hits
312 
313  if( !matchedHit )
314  nonMuonClust.RemoveHit( hit );
315 
316  }// end loop over origSlice hits
317 
318  // Now we have the muon removed original cells ,
319  // Lets calculate the efficiency we're after.
320  int nNonMuonHits = nonMuonClust.NCell();
321  int hitCount = 0;
322 
323  art::PtrVector<rb::CellHit> nonMuonHits = nonMuonClust.AllCells();
324  art::PtrVector<rb::CellHit> mrccSliceHits = mrccSlice->AllCells();
325  rb::SortByPlane(nonMuonHits);
326  rb::SortByPlane(mrccSliceHits);
327 
328  int start = 0;
329  for(int oh = 0; oh < nNonMuonHits; oh++){
330  const art::Ptr<rb::CellHit> hit = nonMuonHits[oh];
331  int plane = hit->Plane();
332  int cell = hit->Cell();
333  int tdc = hit->TDC();
334 
335  for(int h = start; h < nMrCells; h++){
336  const art::Ptr<rb::CellHit> mrHit = mrccSliceHits[h];
337  if( mrHit->Plane() > plane ){
338  start = h;
339  break;
340  }
341  if( mrHit->Plane() == plane &&
342  mrHit->Cell() == cell &&
343  mrHit->TDC() == tdc ){
344  hitCount++;
345  break;
346  }
347  }// end loop over mrcc hits
348  }// end loop over non-muon hits
349 
350  maxEff = (float)hitCount/nNonMuonHits;
351 
352  return matchedSliceIndex;
353 
354  }// End of MatchToOrigSlice
355 
356 
357 
358  /*********************************************************************************************/
359  ////////////////////////////////////////////////////////////////////////
360 
362 
363 }// End of namespace murem
back track the reconstruction to the simulation
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.
Atom< std::string > mrClusterLabel
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
FindParentParams fParams
void produce(art::Event &evt)
unsigned short Plane() const
Definition: CellHit.h:39
A collection of associated CellHits.
Definition: Cluster.h:47
Atom< std::string > clusterLabel
void RemoveHit(const art::Ptr< rb::CellHit > hit)
Remove hit from current cluster.
Definition: Cluster.cxx:290
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
void SetParentIndex(int parentIdx)
Set the index of the parent slice.
Definition: MRCCParent.h:33
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
FindParent(const Parameters &params)
Atom< std::string > mrHitLabel
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
data_t::iterator iterator
Definition: PtrVector.h:60
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
void SetParentEff(float parentEff)
Set the efficiency of MRCC slice wrt the parent.
Definition: MRCCParent.h:39
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: event.h:1
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Atom< std::string > mrccLabel
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
Atom< std::string > hitLabel
void SetParentPur(float parentPur)
Set the purity of MRCC slice wrt the parent.
Definition: MRCCParent.h:45
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
Definition: CellHit.cxx:124
int MatchToOrigSlice(art::Ptr< rb::Cluster > mrSlice, art::Handle< std::vector< rb::Cluster > > origSlices, std::vector< art::Ptr< rb::CellHit > > mrHits, float &maxPur, float &maxEff)