ShowerHitSeparator_module.cc
Go to the documentation of this file.
1 #ifndef SHWHITSEP__H
2 #define SHWHITSEP__H
3 
4 ////////////////////////////////////////////////////////////////////////
5 //
6 // Shower Hit Separator
7 //
8 // leigh.howard.whitehead@cern.ch June 2017
9 //
10 // This is a very simple module that takes the output from the EM
11 // separating CNN and applies a cut on the track-like nature
12 // of each hit. Two hit selections are produced, one containing
13 // track-like hits and one containing shower-like hits. This should
14 // work for any MVA method that associates some sort of weight to
15 // hit objects.
16 //
17 // Can optionally produce an output tree containing the information
18 // to visualise the separation of hits.
19 //
20 ////////////////////////////////////////////////////////////////////////
21 
22 #include <string>
23 #include <vector>
24 #include <utility>
25 #include <memory>
26 #include <iostream>
27 #include <map>
28 
29 // Framework includes
33 #include "art_root_io/TFileService.h"
34 
35 // LArSoft Includes
42 
43 #include "TTree.h"
44 
45 namespace shs {
47  public:
48  explicit ShowerHitSeparator(fhicl::ParameterSet const& pset);
49  void reconfigure(fhicl::ParameterSet const& p);
50  void produce(art::Event& evt);
51  void beginJob(){};
52  void endJob(){};
53 
54  private:
55  template<size_t N> void ReadMVA(std::vector<float> &weights, art::Event& evt);
56  void SaveTree(std::vector<recob::Hit> const& shwHits, std::vector<recob::Hit> const& trkHits);
57 
60  double fMVAOutputCut;
61  bool fSaveTree;
62 
63  std::string fMVAClusterLabel; // Get this from the MVA itself.
64 
65  // Output tree and its variables
66  TTree *fOutTree;
67  Int_t fOutEvent;
68  Int_t fOutRun;
69  Int_t fOutSubrun;
70  Int_t fOutTPC;
71  Int_t fOutCryo;
72  Int_t fOutPlane;
73  Int_t fOutTime;
74  Int_t fOutWire;
75  bool fOutIsShw;
76  bool fOutIsTrk;
77 
79  };
80 
81  // implementation
82 
84 
85  this->reconfigure(pset);
86 
87  // Let HitCollectionCreator declare that we are going to produce
88  // hits and associations with wires and raw digits
91 
92  // Define some histograms if we want the plots
93  if(fSaveTree){
94 
96 
97  fOutTree = tfs->make<TTree>("cnnHitTree","cnnHitTree");
98 
99  fOutTree->Branch("Event",&fOutEvent,"Event/I");
100  fOutTree->Branch("Run",&fOutRun,"Run/I");
101  fOutTree->Branch("Subrun",&fOutSubrun,"Subrun/I");
102  fOutTree->Branch("TPC",&fOutTPC,"TPC/I");
103  fOutTree->Branch("Cryostat",&fOutCryo,"Cryostat/I");
104  fOutTree->Branch("Plane",&fOutPlane,"Plane/I");
105  fOutTree->Branch("Time",&fOutTime,"Time/I");
106  fOutTree->Branch("Wire",&fOutWire,"Wire/I");
107  fOutTree->Branch("IsShower",&fOutIsShw,"IsShower/O");
108  fOutTree->Branch("IsTrack",&fOutIsTrk,"IsTrack/O");
109  }
110  }
111 
113 
114  // Get the parameters from the fcl file
115  fMVALabel = p.get<std::string>("MVALabel");
116  fHitLabel = p.get<std::string>("HitLabel");
117  fMVAOutputCut = p.get<double>("MVAOutputCut");
118  fSaveTree = p.get<bool>("SaveTree");
119 
120  }
121 
122  // Access the information from the MVA method. Very similar to the code in PMAlgTrackMaker_module
123  template<size_t N> void ShowerHitSeparator::ReadMVA(std::vector<float> &weights, art::Event& evt){
124  // Access the MVA. It could have 4 or 3 outputs so try to see which.
126  if (!cluResults){
127  // We don't seem to have an MVA.
128  return;
129  }
130 
131  // This is the best way to get hold of the MVA clusters and tag used to associate the clusters
132  // back to the underlying hit objects.
133  fMVAClusters = cluResults->dataHandle();
134  fMVAClusterLabel = cluResults->dataTag();
135 
136  // Now we have the MVA, we need to access the outputs
137  int trkLikeIdx = cluResults->getIndex("track");
138  int emLikeIdx = cluResults->getIndex("em");
139  if ((trkLikeIdx < 0) || (emLikeIdx < 0)) {
140  // If the variables we want don't exist then we have a problem.
141  return;
142  }
143 
144  // Actually extract the weights
145  const auto & cnnOuts = cluResults->outputs();
146  for (size_t i = 0; i < cnnOuts.size(); ++i){
147 
148  double trkOrEm = cnnOuts[i][trkLikeIdx] + cnnOuts[i][emLikeIdx];
149  double val = 0;
150  if (trkOrEm > 0){
151  // Make sure output is normalised to fall between 0 and 1.
152  val = cnnOuts[i][trkLikeIdx] / trkOrEm;
153  }
154  weights.push_back(val);
155  }
156 
157  }
158 
160 
161  // We want to produce two hit collections
162  recob::HitCollectionCreator showerHits(evt,"showerhits",true,true);
163  recob::HitCollectionCreator trackHits(evt,"trackhits",true,true);
164 
165  // These are the MVA weights
166  std::vector<float> mvaWeights;
167 
168  // Access the MVA. It could have 4 or 3 outputs so try to see which.
169  ReadMVA<4>(mvaWeights,evt);
170  if(mvaWeights.size() == 0){
171  // If we didn't get any weights then try the version with three outputs instead
172  ReadMVA<3>(mvaWeights,evt);
173  }
174 
175  // If we get some weights then separate the hits into two collections
176  if(mvaWeights.size() != 0){
177 
178  const art::FindManyP<recob::Hit> hitsFromClusters(fMVAClusters, evt,fMVAClusterLabel);
179 
180  // We also need the hit selection that was used by the CNN
181  auto inputHits = evt.getHandle<std::vector<recob::Hit> >(fHitLabel);
182  art::FindOneP<raw::RawDigit> rawDigits(inputHits,evt,fHitLabel);
183  art::FindOneP<recob::Wire> recoWires(inputHits,evt,fHitLabel);
184 
185  // At this stage we have the hits from the clusters and the tag for each cluster.
186  // Loop over the clusters and add the hits to the shower-like hit collection.
187  for (size_t c = 0; c != fMVAClusters->size(); ++c){
188 
189  // Get the hits from this cluster
190  auto const& hits = hitsFromClusters.at(c);
191 
192  // Now we must copy each hit and reassociate it to the raw digits and wires
193  for (auto const & h : hits){
194  // The art pointer key gives us the position of "h" in the hit collection
195  recob::Hit thisHit = (*h);
196  auto digit = rawDigits.at(h.key());
197  auto wire = recoWires.at(h.key());
198  // Add the hit and its associated raw digit and wire to the collection.
199  if(mvaWeights[c] < fMVAOutputCut){
200  showerHits.emplace_back(thisHit,wire,digit);
201  }
202  else{
203  trackHits.emplace_back(thisHit,wire,digit);
204  }
205  } // End loop over hits
206  } // End for loop over clusters
207  }
208 
209  std::cout << "CNN output splitter: " << std::endl;
210  std::cout << " - Found " << showerHits.size() << " shower-like hits" << std::endl;
211  std::cout << " - Found " << trackHits.size() << " track-like hits" << std::endl;
212 
213  // Fill the output tree if requested to do so.
214  if(fSaveTree){
215  fOutEvent = evt.event();
216  fOutRun = evt.run();
217  fOutSubrun = evt.subRun();
218  SaveTree(showerHits.peek(),trackHits.peek());
219  }
220 
221  // Put the hit collections and associations into the event
222  showerHits.put_into(evt);
223  trackHits.put_into(evt);
224 
225  }
226 
227  void ShowerHitSeparator::SaveTree(std::vector<recob::Hit> const& shwHits, std::vector<recob::Hit> const& trkHits){
228 
229  // Event, Run and Subrun set outside of this function. Take care of everything else here.
230 
231  for(auto const shwHit : shwHits){
232  fOutTPC = shwHit.WireID().TPC;
233  fOutCryo = shwHit.WireID().Cryostat;
234  fOutPlane = shwHit.WireID().planeID().Plane;
235  fOutTime = shwHit.StartTick();
236  fOutWire = shwHit.WireID().Wire;
237  fOutIsShw = true;
238  fOutIsTrk = false;
239 
240  fOutTree->Fill();
241  }
242 
243  for(auto const trkHit : trkHits){
244  fOutTPC = trkHit.WireID().TPC;
245  fOutCryo = trkHit.WireID().Cryostat;
246  fOutPlane = trkHit.WireID().planeID().Plane;
247  fOutTime = trkHit.StartTick();
248  fOutWire = trkHit.WireID().Wire;
249  fOutIsShw = false;
250  fOutIsTrk = true;
251 
252  fOutTree->Fill();
253  }
254 
255  }
256 
258 
259 } // end of shs namespace
260 #endif
EventNumber_t event() const
Definition: DataViewImpl.cc:85
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
static void declare_products(art::ProducesCollector &collector, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
Definition: HitCreator.cxx:248
ShowerHitSeparator(fhicl::ParameterSet const &pset)
Helper functions to create a hit.
A class handling a collection of hits and its associations.
Definition: HitCreator.h:508
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
void SaveTree(std::vector< recob::Hit > const &shwHits, std::vector< recob::Hit > const &trkHits)
T get(std::string const &key) const
Definition: ParameterSet.h:271
void ReadMVA(std::vector< float > &weights, art::Event &evt)
p
Definition: test.py:223
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
void emplace_back(recob::Hit &&hit, art::Ptr< recob::Wire > const &wire=art::Ptr< recob::Wire >(), art::Ptr< raw::RawDigit > const &digits=art::Ptr< raw::RawDigit >())
Adds the specified hit to the data collection.
Definition: HitCreator.cxx:290
RunNumber_t run() const
Definition: DataViewImpl.cc:71
void put_into(art::Event &)
Moves the data into an event.
Definition: HitCreator.h:652
ProducesCollector & producesCollector() noexcept
size_t size() const
Returns the number of hits currently in the collection.
Definition: HitCreator.h:636
Declaration of signal hit object.
void reconfigure(fhicl::ParameterSet const &p)
art::Handle< std::vector< recob::Cluster > > fMVAClusters
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
TCEvent evt
Definition: DataStructs.cxx:7
std::vector< recob::Hit > const & peek() const
Returns a read-only reference to the current list of hits.
Definition: HitCreator.h:665
static std::unique_ptr< MVAReader > create(const art::Event &evt, const art::InputTag &tag)
Definition: MVAReader.h:110
QTextStream & endl(QTextStream &s)