FindMREParent_module.cc
Go to the documentation of this file.
1 // C/C++ & root includes
2 #include <iostream>
3 #include "vector"
4 #include <map>
5 
6 // ART includes
7 
21 
22 //NOvA includes
23 #include "Utilities/AssociationUtil.h"
26 #include "ReMId/ReMId.h"
27 #include "RecoBase/Cluster.h"
28 #include "MCCheater/BackTracker.h"
29 
30 namespace murem
31 {
32 
33  namespace
34  {
35  struct FindMREParentParams
36  {
37  template<class T> using Atom = fhicl::Atom<T>;
38  template<class T> using Table = fhicl::Table<T>;
39  using Comment = fhicl::Comment;
40  using Name = fhicl::Name;
41 
42  Atom<std::string> mrccLabel
43  {
44  Name("MRCCLabel"),
45  Comment("Label of the mrcc raw digit module")
46  };
47  Atom<std::string> mreLabel
48  {
49  Name("MRELabel"),
50  Comment("Label of the mre raw digit module")
51  };
52  Atom<std::string> mreClusterLabel
53  {
54  Name("MREClusterLabel"),
55  Comment("Label of mre slicer module")
56  };
57  Atom<std::string> clusterLabel
58  {
59  Name("ClusterLabel"),
60  Comment("Label of the parent slicer module")
61  };
62 
63  };// end of parameters' struct
64 
65  }
66 
68  {
69  public:
70  // Allows 'nova --print-description' to work
72 
73  explicit FindMREParent(const Parameters& params);
74  virtual ~FindMREParent();
75 
76  void produce(art::Event& evt);
77 
78  private:
79  FindMREParentParams fParams;
80 
81  };
82 
83  /*************************************************************/
84  ///////////////////////////////////////////////////////////////
85 
87  : fParams(params())
88  {
89 
90  produces<std::vector<murem::MRCCParent> >();
91  produces< art::Assns<murem::MRCCParent, rb::Cluster> >();
92  }
93 
94 
95  /*************************************************************/
96  ///////////////////////////////////////////////////////////////
97 
98 
100  {
101  }
102 
103  /*************************************************************/
104  ///////////////////////////////////////////////////////////////
105 
107  {
108  // the mreparent object to be put in event
109  std::unique_ptr< std::vector< murem::MRCCParent > >
110  mrccParents(new std::vector< murem::MRCCParent >);
111 
112  // the mreparent object to slice association
113  std::unique_ptr< art::Assns< murem::MRCCParent,rb::Cluster > >
115 
116  std::vector< murem::MRCCParent> mrccParentsVec;
117 
118  // mre slices
119  art::Handle< std::vector<rb::Cluster> > mreslicehandle;
120  evt.getByLabel(fParams.mreClusterLabel(), mreslicehandle);
121  art::PtrVector<rb::Cluster> mreslices;
122  for(unsigned i=0; i < mreslicehandle->size(); i++){
123  art::Ptr<rb::Cluster> ptr(mreslicehandle, i);
124  mreslices.push_back(ptr);
125  }
126 
127  // mre electron mc truth
129  evt.getByLabel(fParams.mreLabel(), elechandle);
130  std::vector<art::Ptr<simb::MCTruth> > elecs;
131  art::fill_ptr_vector(elecs, elechandle);
132 
133  // removed muon particles
135  evt.getByLabel(fParams.mrccLabel(), muonhandle);
136  std::vector<art::Ptr<simb::MCParticle> > muons;
137  art::fill_ptr_vector(muons, muonhandle);
138 
139  // parent slices
141  evt.getByLabel(fParams.clusterLabel(), slicehandle);
142  std::vector<art::Ptr<rb::Cluster> > slices;
143  art::fill_ptr_vector(slices, slicehandle);
144 
146  std::vector<std::vector<cheat::NeutrinoEffPur> > truthtable =
147  bt->SlicesToMCTruthsTable( mreslices);
148 
149  if(!truthtable.empty()){
150 
151  // Map the mre slice to an electron mctruth
152  // electron mctruths are associated with
153  // an mcparticle that corresponds to the
154  // removed muon, which is asociated with
155  // the parent event slice.
156 
157  for(unsigned int isli = 0; isli < mreslices.size(); isli++){
158 
159  if(mreslices[isli]->IsNoise())
160  continue;
161 
162  if(truthtable[isli].empty())
163  continue;
164 
165  // Find the best match for this slice
166  int bestidx = -1;
167  float besteff = 0;
168  for(size_t it = 0; it < truthtable[isli].size(); ++it){
169  if( truthtable[isli][it].efficiency > besteff){
170  bestidx = it;
171  besteff = truthtable[isli][it].efficiency;
172  }
173  }
174 
175  if( bestidx < 0 ||
176  truthtable[isli].at(bestidx).efficiency <= 0 )
177  continue;
178 
179  art::Ptr< simb::MCTruth > electruth =
180  truthtable[isli].at(bestidx).neutrinoInt;
181 
183  fop({electruth}, evt, fParams.mreLabel());
184 
185  if(!fop.isValid() )
186  continue;
187 
188 
189  // the electron truth is associated with
190  // the removed muon mc particle. Grab that
191  art::Ptr<simb::MCParticle> muon = fop.at(0);
192 
194  fos({muon}, evt, fParams.mrccLabel());
195 
196  if(!fos.isValid() && fos.at(0))
197  continue;
198 
199  art::Ptr<rb::Cluster> parentslice = fos.at(0);
200  int sliceIdx = -5;
201  for(unsigned int i = 0; i < slices.size(); i++){
202  if(slices[i] == parentslice)
203  sliceIdx = i;
204  }
205 
206  if(sliceIdx == -5 )
207  continue;
208 
209  // now find the efficiency and purity numbers
210  MRCCParent mrccParent;
211  mrccParent.SetParentPur( 1 );
212  mrccParent.SetParentEff( 1 );
213  mrccParent.SetParentIndex( sliceIdx );
214 
215  mrccParents->push_back( mrccParent );
216  util::CreateAssn( *this, evt, *mrccParents,
217  mreslices[isli] ,*assns);
218 
219  }// end loop over mre slices
220  }// mre truth table is not empty
221 
222  evt.put( std::move(mrccParents) );
223  evt.put( std::move(assns) );
224  }// End of produce method
225 
226 
228 
229 }// End of namespace murem
FindMREParent(const Parameters &params)
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.
set< int >::iterator it
Atom< std::string > mreClusterLabel
Atom< std::string > clusterLabel
DEFINE_ART_MODULE(TestTMapFile)
Atom< std::string > mreLabel
Particle class.
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
void SetParentIndex(int parentIdx)
Set the index of the parent slice.
Definition: MRCCParent.h:33
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void produce(art::Event &evt)
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
size_type size() const
Definition: PtrVector.h:308
FindMREParentParams fParams
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
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
void SetParentPur(float parentPur)
Set the purity of MRCC slice wrt the parent.
Definition: MRCCParent.h:45
void efficiency()
Definition: efficiency.C:58
Atom< std::string > mrccLabel